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respect  to  existing  oxidation  mechanism  reduction  schemes.  The  manuscripts  are  individually 
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EXECUTIVE  SUMMARY 


The  highlights  of  the  results  from  the  manuscripts  in  Appendices  1-4  are  here  summarized. 

The  challenge  of  modeling  turbulent  reactive  flows  is  so  considerable  that  the  activity  has 
traditionally  been  decomposed  into  its  two  essential  parts:  kinetics  and  turbulence.  Usually, 
modeling  of  chemical  kinetics  has  proceeded  on  a  separate  path  from  that  of  turbulence  which 
also  includes  canonical  models  for  turbulence/reaction  interaction.  The  only  constraint  to  kinetic 
modeling  was  that  it  should  be  compact  enough  to  be  computationally  efficient  when  included  in 
a  complex  turbulent  combustion  code.  However,  there  are  definite  advantages  on  approaching 
chemical  kinetic  modeling  in  a  similar  manner  to  turbulent  flow  modeling  because  if  the 
concepts  are  similar,  the  hope  is  that  the  models  will  mesh  better  and  the  results  will  be  easier  to 
understand.  This  was  the  approach  taken  in  this  study.  The  spirit  of  the  chemical  kinetic 
modeling  approach  is  that  of  Large  Eddy  Simulations  (LES)  in  which  kinematically 
energetically-signi (leant  flow  scales  are  computed  and  the  others  are  modeled;  in  turbulence,  the 
large  flow  scales  constitute  the  fonner  category  and  the  small  flow  scales  constitute  the  latter 
category.  The  chemical  kinetics  parallel  is  to  obtain  a  model  in  which  one  retains  only  the 
thermodynamically  energetically-significant  chemical  scales  as  progress  variables,  and  models 
the  fate  of  the  other  scales.  But  the  parallel  approach  between  kinetics  and  turbulence  was 
extended  even  further.  In  turbulence  modeling,  a  common  methodology  is  to  assess  the  behavior 
of  the  modeled  scales  by  analyzing  databases  created  using  Direct  Numerical  Simulations  (DNS) 
in  which  all  flow  scales  are  computed;  indeed,  current  experimental  diagnostics  do  not  pennit 
the  same  thoroughness  of  infonnation  as  that  obtained  from  DNS.  The  DNS  data  is  analyzed  in 
what  is  called  an  a  priori  study  to  inquire  about  the  behavior  of  the  small  scales  and  propose 
mathematical  fonns  which  fit  this  behavior.  The  a  priori  study  is  followed  by  an  a  posteriori 
study  where  the  proposed  mathematical  forms  are  inserted  into  the  model  to  evaluate  its 
performance  when  compared  to  the  DNS  database  at  the  LES  resolution. 

A  complete  analogy  between  turbulence  and  kinetics  was  made  by  observing  that  kinetic 
elemental  or  skeletal  mechanisms  can  serve  for  reduced  kinetics  the  role  that  DNS  serves  for 
LES,  in  which  case  reduced  kinetic  mechanisms  can  be  viewed  as  the  complement  to  LES  in 
achieving  the  goal  of  accurate  computationally-efficient  turbulent  reactive  flow  simulations.  The 
analogy  between  reduced  kinetic  models  and  LES  is  not  entirely  surprising  since  each  chemical 
species  has  a  characteristic  time  scale  and  in  the  kinetic  reduction  it  is  desirable  to  compute  only 
those  entities  (e.g.  species,  combination  of  species,  radicals,  combination  of  radicals,  etc.)  having 
essential  characteristic  time  scales  (to  be  defined)  and  model  the  kinetics  of  the  remaining 
entities. 

Thus,  there  were  two  important  components  to  this  study,  namely  the  a  priori  analysis  and 
the  a  posteriori  evaluation.  Whereas  turbulence  modeling  benefits  from  decades  of  work  using 
the  DNS/LES  concept  which  originated  in  atmospheric  turbulence  predictions  in  the  1960s,  the 
present  work  is  the  first  investigation  to  take  this  approach  in  kinetic  reduction  modeling. 
Therefore,  it  was  first  necessary  to  produce  a  categorization  of  scales  analogous  to  the  large  and 
small  scales  of  turbulence,  then  it  was  required  to  propose  mathematical  forms  for  the  scales  that 
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will  be  modeled  rather  than  computed  as  progress  variables  in  the  reduced  kinetic  model,  and 
finally  it  was  required  to  perform  an  a  posteriori  study  in  order  to  evaluate  the  chemical  kinetic 
model  versus  the  elemental  or  skeletal  mechanism  for  those  species  predicted  by  the  reduced 
model.  The  present  model  depicts  a  constant-volume  situation,  so  as  to  be  consistent  with  the 
requirement  of  a  LES  grid. 

We  have  proposed  such  a  categorization  [1-4]  through  the  definition  of  a  total  constituent 
molar  density  which  is  a  progress  variable  representing  the  heavy  species  kinetics,  and  through 
the  partition  of  the  light  species  set  into  a  set  of  modeled  quasi-steady  species  and  a  set  of 
progress  variable  species.  By  definition,  constituent  radicals  are  those  composing  species  with  a 
carbon  number  larger  than  or  equal  to  3.  These  species  are  called  ‘heavy’,  and  their  complement 
in  the  ensemble  of  species  is  called  Tight’.  The  total  constituent  molar  density  is  the  sum  of  the 
individual  constituents’  molar  densities.  Thus,  rather  than  following  all  species  through  their 
reaction  coordinates,  we  follow  a  reduced  set  of  reaction  coordinates  (i.e.  progress  variables); 
this  reduced  set  is  called  a  base.  The  a  priori  model,  with  numerous  results  was  described  in  [1- 
4].  The  behavior  of  the  base  set  was  examined  a  priori  in  the  context  of  a  detailed  reaction  model 
for  n-heptane  combustion  as  given  by  LLNL  [5],  and  the  proposed  mathematical  forms  were 
tested  against  the  LLNL  kinetic  scheme.  The  first  a  posteriori  evaluation  of  the  model  was 
presented  in  [4]  with  the  goal  of  inquiring  whether  the  qualitative  aspects  of  the  reaction 
evolution  were  captured  and  identifying  refinements  that  could  be  made  to  improve,  if  necessary, 
the  quantitative  agreement  with  the  LLNL  database. 

To  give  perspective  to  our  model,  one  must  consider  that  there  are  several  ways  of 
simplifying  the  chemistry,  but  the  tenninology  employed  for  qualifying  this  simplification  is  not 
unified  and  instead  depends  on  the  author.  For  example,  in  [6]  the  ‘complete  mechanism’  is 
denoted  as  that  in  which  all  kinetic  processes  have  been  taken  into  account;  a  ‘detailed’  or 
‘skeletal’  mechanism  is  that  where  the  complexity  of  the  complete  mechanism  has  been  reduced 
for  a  given  purpose,  while  keeping  all  main  kinetic  pathways  important  for  the  situation  at  hand 
(reactions  considered  unimportant  through  sensitivity  analysis  are  eliminated  and  possibly 
lumping  of  isomers  decreases  the  number  of  species,  and  thus  decreases  the  number  of 
differential  equations);  a  ‘semi-global’  mechanism  is  defined  as  that  involving  5-10  species  and 
in  which  most  chemical  pathways  have  been  neglected;  and  a  ‘single-step’  mechanism  is  that 
where  all  intennediary  species  have  been  removed  and  only  reactants  and  products  appear.  The 
authors  of  [6]  acknowledge  that  the  boundaries  between  defined  categories  are  not  rigid. 
According  to  their  definition,  the  skeletal  mechanism  represents  the  lower  limit  of  the  detailed 
mechanism.  Another  frequent  terminology  in  kinetic  modeling  is  that  of  a  ‘reduced’  mechanism, 
meaning  that  the  mechanism  is  smaller  than  the  original  one.  For  example,  there  is  a  large 
number  of  smaller-than-the-LLNL-mechanism  for  n-heptane,  from  that  of  [7]  that  is  a  30-step 
mechanism,  to  that  of  [8]  having  a  set  of  14  reactions  and  16  species,  to  that  further  adjusted  in 
[9]  to  a  set  of  4  steps  and  5  species.  Recently,  researchers  have  announced  [10]  a  new  n-heptane 
skeletal  mechanism  of  78  species  with  359  or  317  reactions,  another  skeletal  mechanism  with  68 
species  and  283  reactions,  and  a  reduced  mechanism  with  55  species  and  51  global  steps  (i.e. 
which  can  be  selected  from  either  elementary  or  lumped  reactions,  their  rate  expressions  being 
the  linear  combination  of  the  rates  of  the  elementary  reactions.).  This  means  that  the  minimum 
number  of  progress  variables,  that  is,  the  number  of  variables  for  which  differential  equations  are 
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solved,  is  55.  A  remaining  number  of  13  species  is  determined  using  algebraic  equations  by 
invoking  the  quasi-steady  state  assumption.  All  these  mechanisms  were  tested  in  the  range  of 
latm  to  50atm  for  the  pressure,  600K  to  1800K  for  the  ignition  temperature,  and  0.5  to  1.5  for 
the  equivalence  ratio;  testing  involved  comparisons  with  the  full  kinetic  mechanism  [5],  and  also 
with  flame  data. 

There  are  several  kinetic  mechanism  methodologies  which  can  be  generally  classified  as 
either  ‘fixed’  or  ‘on-the-fly’.  In  fixed  methodologies,  the  kinetic  mechanism  is  smaller  than  the 
original  one  and  when  entering  a  reaction  calculation  the  same  mechanism  is  meant  to  be  used 
during  the  entire  calculation.  An  example  of  a  fixed  mechanism  is  given  in  [10].  In  contrast, 
for  ‘on-the-fly’  mechanisms,  criteria  are  used  to  sample  during  a  reaction  calculation  the  entire 
mechanism,  and  at  each  time  a  reduction  is  found,  which  may  vary  during  the  calculation.  An 
example  of  such  a  method  is  that  of  the  intrinsic  low-dimensional  manifold  (ILDM)  [11].  The 
ILDM  method  is  based  on  the  observation  that  in  a  chemically  reactive  mixture  at  constant 
pressure,  constant  enthalpy  and  constant  element  composition,  the  path  of  reaction  in  the  high¬ 
dimensional  state  space  lies  in  low-dimensional  manifolds  far  from  the  final  equilibrium  point; 
these  manifolds  are  intrinsic,  meaning  that  except  for  the  degree  of  reduction,  the  reaction 
mechanism  in  tenns  of  Arrhenius-type  elementary  reactions  is  the  only  source  of  infonnation 
necessary  to  build  the  approximation.  It  was  found  [11]  that  ILDM  is  not  easily  implementable 
in  some  situations.  Particularly,  it  (i)  has  the  drawback  of  becoming  very  difficult  to  handle 
with  increasing  number  of  C  atoms  in  the  fuel  [11],  and  we  note  that  heavy  species  (i.e.  large 
number  of  C  atoms)  enter  the  composition  of  most  practical  fuels,  and  (ii)  does  not  work  at  low 
temperatures  (the  ignition  region),  or  very  lean  or  rich  mixture  conditions,  which  are  precisely 
the  conditions  where  a  very  large  number  of  pollutants  are  formed. 

In  our  methodology  [1-4],  a  much  larger  reduction  than  that  obtained  in  [10]  has  been 
achieved  because  instead  of  the  55  progress  variables  based  on  species,  as  in  [10],  we  only  have 
only  12  progress  variables  based  on  constituent  radicals  and  light  species,  as  defined  below. 
Moreover,  instead  of  the  n-heptane  mechanism  being  only  valid  in  the  range  of  latm  to  50atm, 
600K  to  1800K  and  0.5  to  1.5  for  the  equivalence  ratio,  our  range  is  5bar  to  60bar,  temperatures 
as  low  as  600K  and  as  high  as  2500K  depending  on  the  pressure,  and  0.125  to  8  range  for  the 
equivalence  ratio.  The  slightly  higher  pressure  is  helpful  for  diesel,  gas  turbine  and  HCCI 
engines;  the  higher  temperature  at  the  upper  limit  is  important  for  NOx  production  reactions  to  be 
grafted  on  the  present  mechanism  in  future  work;  the  considerably  wider  equivalence  ratio 
regime  is  crucial  for  NOx  production  in  the  lean  regime  and  CO  and  soot  production  in  the  rich 
regime  (future  work). 

Moreover,  our  model  is  hierarchical  by  construct,  meaning  that  it  is  naturally  extendable  to 
higher  carbon-number  hydrocarbons  (a  primary  necessity  for  modeling  iso-octane  and  ring 
hydrocarbons),  thus  having  advantages  in  this  respect  over  the  capabilities  of  ILDM.  The 
hierarchical  aspect  comes  in  because  the  progress  variables  are  not  necessarily  species,  but  also 
include  species  ‘constituents’  [1-4],  defined  very  much  like  in  group  additivity  theory  [12].  In 
contrast  to  groups  [12],  constituents  take  into  account  only  interactions  with  adjacent  groups  and 
only  first  order  (compositional)  effects.  Examination  of  the  LLNL  n-heptane  mechanism  over  a 
wide  range  of  pressures  and  stoichiometries  showed  that  most  of  the  constituent  rates  are  quasi- 
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steady;  exceptions  are  only  for  radicals  with  minor  contribution  to  the  total.  Our  modeling  effort 
reduced  the  LLNL  n-heptane  mechanism  from  160  species  (progress  variables)  and  1540 
reactions  to  12  progress  variables  (11  light  species  and  the  global  constituent  mole  fraction)  for 
which  16  rates  (the  contributions  of  the  heavy  species  to  the  light  species  -11  gain  rates  and  4 
loss  rates-  and  the  global  constituent  rate)  have  been  modeled  by  fitting,  162  conventional 
reaction  rates  (light  species  interacting  among  themselves)  and  1 1  other  functional  forms  (i.e.  fits 
for  the  heat  capacity  at  constant  pressure,  the  enthalpy  release  rate  of  the  heavy  species,  and  the 
molar  fraction  of  quasi-steady  light  species).  The  reduced  model  comprises  16  rate  functions,  1 1 
curves  for  light  quasi-steady  mole  fractions  and  the  heat  capacities  of  the  heavies  fitted  in  terms 
of  a  normalized  temperature,  the  stoichiometric  ratio,  the  initial  temperature,  and  an  initial  scaled 
pressure  parameter.  This  is  a  reduction  of  the  number  of  progress  variables  by  an  order  of 
magnitude  from  the  LLNL  elementary  mechanism.  This  compact  model  was  verified  by 
comparing  temperature  evolution  paths  from  calculations  using  the  LLNL  dataset  to 
corresponding  paths  from  the  reduced  model  of  constituents  and  light  species;  deviations  were  at 
most  a  few  percent.  Therefore,  going  from  n-heptane  to  heavier  hydrocarbons  the  number  of 
constituents  may  increase  (e.g.  with  aromatics),  but  if  as  for  n-heptane  they  can  all  be  bundled 
into  a  single  total  constituent  mole  fraction  that  is  the  single  representative  progress  variable  of 
the  constituents,  it  is  predicted  that  the  computational  time  should  remain  essentially  the  same. 

The  detailed  analysis  of  the  database  had  an  additional  advantage:  the  finding  of  a  scaling  of 
the  global  constituents  mole  fraction  and  of  a  temperature  non-dimensional  variable  that  allowed 
the  establishment  of  a  self-similarity  of  the  global  constituents’  mole  fraction  with  respect  to  the 
non-dimensionalized  temperature,  pressure,  equivalence  ratio  and  temperature  at  which  the 
global  constituents  reaction  rate  becomes  positive  [4].  This  self  similarity  is  a  reduction  in 
problem  dimensionality  but  is  not  to  be  confused  with  ILDM  which  was  for  species,  whereas  the 
present  one  is  for  the  constituents.  The  fact  that  the  constituents  can  embody  the  evolution  of  the 
heavy  species  is  physically  consistent  with  the  fact  that  as  the  temperature  increases,  the  heavy 
species  decompose  and  no  longer  play  a  role;  instead,  the  products  of  this  decomposition 
determine  the  reaction  evolution. 

Finally,  the  preliminary  a  posteriori  study  revealed  that  the  qualitative  features  of  the 
reaction  evolution  are  well  captured  but  that  fitting  of  the  reaction  rates  in  the  initial  part  of  the 
incubation  region  (where  the  non-dimensional  temperature  is  smaller  than  10'  )  must  be 
improved  to  obtain  good  quantitative  results.  In  this  region,  the  temperature  increases  by 
approximately  5K  (we  call  this  the  quiescent  part  of  the  incubation  region),  before  the  reaction 
transitions  to  the  next  part  of  the  incubation  region  (which  we  call  the  active  part)  where  the 
temperature  increases  to  approximately  1000K  which  results  in  the  alteration  of  the  heavy 
species  structure  and  brings  about  their  decomposition.  This  improved  fitting  has  been 
performed  and  is  currently  being  evaluated  a  posteriori  in  the  renewed  grant  from  ARO  on  this 
subject. 
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A  new  method  in  modeling  and  simulations  of  complex 

oxidation  chemistry 
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Jet  Propulsion  Laboratory** ,  California  Institute  of  Technology, 
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A  simplified  model  is  proposed  for  the  kinetics  of  alkane  oxidation  in  air,  based  on  a 
decomposition  of  heavy  (carbon  number  >  3)  hydrocarbons  into  a  13  constituent  radical 
base.  The  behavior  of  this  base  is  examined  in  test  computations  for  heptane  utilizing 
Chemkin  II  with  LLNL  data  inputs.  Emphasis  is  placed  on  prediction  of  the  heat  release  and 
temperature  evolution.  At  stoichiometric  conditions,  the  total  constituent  molar  density 
was  found  to  follow  a  quasi-steady  rate  which  is  a  simplification  in  the  modeling  of  its 
reaction  rate. 


Introduction 

Computational  models  for  turbulent  reactive  flows  present  challenges  that  have  yet  to  be  fully  met. 
First,  combustion  of  typical  hydrocarbons  proceeds  through  the  formation  of  hundreds  of  species  interacting 
through  thousands  of  reactions.  Second,  a  complete  numerical  resolution  of  high  Reynolds  number  flow 
requires  a  very  large  number  of  grid  points  that  is  beyond  current  typical  supercomputer  capabilities.  Thus, 
combining  detailed  reaction  rate  models  with  the  fundamental  flow  equations  yields  a  computational  problem 
that  is  computationally  prohibitive.  A  drastic  (at  least  order  of  magnitude)  reduction  in  the  dimension  of 
reaction  kinetics  is  required  for  a  feasible  computational  model.  In  spirit,  this  reduction  is  similar  to  that 
to  be  achieved  in  Large  Eddy  Simulations  (where  energetically  significant  flow  scales  are  computed  and  the 
others  are  modeled)  with  respect  to  Direct  Numerical  Simulations  (where  all  flow  scales  are  computed). 
Since  each  species  has  a  characteristic  time  scale,  in  the  kinetic  reduction  it  is  also  desirable  to  compute 
only  those  species  that  have  an  essential  characteristic  time  scale  (in  a  sense  to  be  defined)  and  model  the 
kinetics  of  the  remaining  species. 

We  categorize  all  species  involved  in  hydrocarbon  oxidation  into  light  species  and  heavy  species;  the 
heavy  species  are  those  having  a  carbon  number  n  ^  3.  Based  on  this  definition,  a  method  for  simplifying 
alkane  combustion  is  proposed  that  focusses  on  a  small  set  of  constituent  radicals  (concept  to  be  defined  in 
Section  I)  of  the  heavy  species  rather  than  on  the  much  larger  number  of  these  species.  That  is,  rather  than 
following  all  species  through  their  reaction  coordinates,  we  will  follow  a  reduced  set  of  reaction  coordinates; 
this  reduced  set  is  called  a  base.  Our  primary  objective  is  the  accurate  determination  of  the  energetics, 
i.e.  heat  release  and  temperature  evolution.  In  Section  I,  we  define  the  constituent  radicals  and  describe 
the  formation  of  the  constituent  radical  base.  In  Section  II,  the  behavior  of  the  base  set  is  examined  in 
the  context  of  a  detailed  reaction  model  for  heptane  combustion  as  given  by  LLNL  investigators.1  Finally, 
Section  III  is  devoted  to  a  discussion  on  how  to  use  curve  fits  developed  in  Section  II  that  describe  rate  and 
species  behavior  and  to  a  perspective  on  the  further  work  required  to  complete  a  reduced  model  for  heptane 
(and  other  alkane)  combustion. 
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I.  Representation  of  species  through  a  reduced  base  set 


A.  Heavy  species  decomposition  into  constituent  radicals 

Aside  from  the  molar  reaction  rates,  determining  the  temperature  evolution  in  a  reactive  system  requires 
knowledge  of  the  species  molar  enthalpies  and  heat  capacities.  For  species  i,  the  molar  enthalpy  may  be 
expressed  as 

hi  =  h°i+hi{T)  (1) 

where  h ?  is  the  heat  of  formation,  hi  is  the  sensible  enthalpy  and  T  is  the  temperature.  The  heat  of  formation 
is  taken  at  reference  conditions,  pref  =  1  bar  and  Tref  =  25  C  =  298.15  K  where  p  is  the  pressure,  giving 

hi=  [  Cp'idT  (2) 

J  Tref 

where  CPti  is  the  molar  constant-pressure  heat  capacity.  A  heat  of  combustion  for  species  i  is  given  by 

hi  =hi  ~Ylwiih°j  (3) 

j 

where  the  index  j  denotes  the  species  set  of  air  (02,N2)  and  final  combustion  products  (H20,C02).  Equation 
3  is  based  on  water  remaining  in  vapor  state  and  not  condensing  into  liquid.  Consider  the  count  of  atoms 
(H,C,0,N)  in  species  i.  Hypothetically,  we  can  consider  species  i  as  being  composed  of  the  air  and  final 
combustion  product  species.  By  definition,  Wji  is  the  number  of  species  j  in  species  i  that  produces  the 
correct  atom  count.  The  molar  densities  of  species  in  this  set  are 

Ai  =  Ai  »'./^  (4) 

i¥=3 

where  N®  is  the  effective  total  molar  density  of  species  j.  The  same  exercise  of  atom  counting  and  species 
representation  by  the  j  set  can  be  performed  for  all  species  in  the  mixture  to  obtain  N(- .  That  is,  the  atomic 
elements  molar  densities  are  represented  by  species  j  molar  densities.  Conservation  of  atoms  in  the  reactions 
makes  unnecessary  utilizing  rate  expressions  for  reaction  coordinates  of  air  species  and  complete-combustion 
products,  as  a  one-to-one  map  exists  between  these  and  the  atom  set. 

Using  information  from  the  NIST  Chemistry  WebBook  website,2  CRC  Handbook  of  Chemistry  and 
Physics,3  the  GRI  website4  and  the  NASA  Glenn  website,5  plots  of  h\  for  alkanes  and  also  for  alkenes 
having  the  carbon  double  bond  at  the  molecular  chain  end  (i.e.  vinyl  radical)  show  a  linear  variation  with  n 
(not  shown).  Further,  at  a  fixed  value  of  T,  the  corresponding  Cpp  also  vary  linearly  with  n.  Both  h\  and  Cpp 
linear  variation  with  n  is  functionally  very  accurate.  This  implies  that  these  thermodynamic  properties  may 
be  considered  as  determined  by  given  properties  of  constituent  radicals  CH2,  CH3  and  C2H3  that  together 
form  these  hydrocarbons.  That  is,  each  heavy  species  is  decomposed  into  constituent  radicals  which  are 
here  used  as  their  building  blocks.  This  decomposition  is  consistent  with  the  group  additivity  concept  of 
Benson0  (see  also  Reid  et  al.' );  however,  here  it  is  desired  to  have  forms  with  greater  simplicity  in  which  only 
first  order  (compositional)  effects  are  allowed.  For  hydrocarbon  species  with  double  bonds  not  at  the  chain 
end  or  for  oxidized  species,  the  simplicity  of  constituent  radicals  as  a  base  no  longer  rigorously,  but  only 
approximately,  holds.  To  calculate  h\  once  either  breakup  of  a  hydrocarbon  into  free  radicals  or  oxidation 
reactions  occur,  the  assumption  is  made  that  one  can  associate  a  fixed  partial  enthalpy  with  each  constituent 
radical,  independent  of  the  species  to  which  they  may  belong.  The  information  from  the  cited  references 
shows  that  although  this  assumption  is  not  entirely  correct  (when  subtracting  contributions  from  CH2,  CH3 
and  C2H3),  the  variation  across  species  is  not  substantial;  as  a  result,  an  average  over  species  is  here  taken. 
The  error  in  the  resulting  h\  is  rather  small  (not  shown);  two-digit  accuracy  is  typically  reached  (with  best 
accuracy  for  the  dominant  species). 

Plots  of  Cp/Ru  versus  logT,  where  Ru  is  the  universal  gas  constant,  are  also  nearly  linear  for  the 
constituent  radicals,  as  shown  in  Fig.  1.  Thus,  for  this  base,  the  form 

5CP/RU  =  ah  +  bh  MT/Tref)  (5) 

is  assumed  with  ah  and  bh  being  constant.  Properties  of  any  constituent  radical  in  a  stable  molecule  having 
valence  bonds  satisfied  are  not  the  same  as  when  the  radical  is  in  free  form,  i.e.  with  one  or  more  available 
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bonds.  In  the  same  manner,  properties  of  heavy  radicals  are  different  from  the  sum  of  their  constituent 
radicals.  For  example,  the  difference  in  heat  of  formation  between  constituent  and  free  radical  is  defined  as 

6hf  =  /i£ee  radical  -  Sh°  (6) 

where  Sh°  is  the  constituent  radical  partial  heat  of  formation.  For  the  set  of  13  constituent  bases  listed  in 
Table  1,  values  of  6h° ,6hc,6h/ ,  ah  and  bh  are  given  in  that  table.  The  heavy  species  that  are  free  radicals 
also  require  5hj  and  6CP  corrections  (i.e.  additional  ah  and  bh  constants)  due  to  unsatisfied  valence  bonds. 
These  corrections  may  be  estimated  using  increments  in  hydrogen  bond  strengths  given  by  Lay  et  al.8  along 
with  information  from  LLNL1  to  NASA  Glenn.5  As  an  example,  let  R  denote  a  radical  produced  by  hydrogen 
abstraction;  RH  is  the  parent  and  Eb  is  the  bond  energy.  Then, 


h°(RH)  +  Eb  =  h°(R)  +  h°(H). 

(7) 

Also, 

h°(RH)  =  Sh°(R)  +  5h°(H) 

(8) 

giving 

6hf(R)  =  Eb  +  6h°(H)  -  h°(H) 

(9) 

with  Sh°(H)  ~  (10?rb  —  33)  kJ/mol  where  nb  is  the  number  of  valence  bonds  of  C  atoms  attached  to  the 
C  atom  of  abstraction  (e.g.  nb  =  0  for  methane,  nb  =  3  for  alkynes).  The  rules  that  were  devised  to 
estimate  Sh-f  and  SCP  for  heavy  radicals  are  given  in  Table  2  along  with  necessary  corrections  for  a  few 
stable  molecules. 

B.  Light  species  as  part  of  the  base 

The  light  species  are  not  subject  to  meaningful  decomposition.  Values  of  hc  from  the  literature,  Cp  s  fitted 
according  to  eq.  5  and  other  properties  of  the  light  molecules  and  light  free-form  of  the  constituent  radicals 
are  given  in  Table  3. 

C.  Final  base  set 

The  reaction  base  set  is  composed  of  13  constituent  radicals  (Table  1)  and  26  light  molecules  or  free  radicals 
(Table  3).  The  entire  base  set  may  not  be  needed  for  a  particular  situation.  Furthermore,  the  atomic  balance 
eliminates  the  need  to  calculate  4  of  the  light  species. 

II.  Base  behavior  in  stirred  reactor  heptane  combustion 

The  evolution  of  T  and  the  molar  densities  N;t  was  determined  in  test  cases  by  utilizing  the  Chemkin 
II  software  package  along  with  species  and  rate  data  from  the  LLNL  website.1  This  data  encompasses  160 
species  interacting  through  1540  reaction  rates  (770  forward  and  770  corresponding  backward).  The  base  set 
excludes  dicarbon  and  also  the  light  species  (from  Table  3)  C,  C2,  N,  NO  and  NO2  as  they  are  not  necessary 
(e.g.  if  the  focus  were  on  soot  or  NO^,  some  of  these  species  would  enter  the  base  set).  Calculations  were 
performed  at  initial  pressures  po  from  1  to  60  bar.  Most  calculations  were  performed  for  a  stoichiometric 
mixture;  a  smaller  number  of  calculations  were  made  for  molar  equivalence  ratios  <j>  of  0.5  and  2.0. 

To  test  the  compatibility  of  using  the  constituent  properties  for  finding  the  T  evolution,  the  energy 
equation 

NC>  f  =  -!>-«*.  R,S{^)reacUm  (10) 

where  N  =  JV  N{  and  NCP  =  NiCPti,  was  solved  twice  and  simultaneously  with  different  enthalpy  inputs. 
In  both  cases,  Ri  is  determined  from  the  full  rate  LLNL1  model;  Cp  and  hi  were  in  one  case  determined 
through  the  LLNL  model  (which  contains  the  enthalpies)  and  in  the  other  case  through  the  simplified  forms 
given  in  Tables  1  through  3.  The  deviations  between  the  two  T  profiles  values  were  typically  no  more  than 
0.3%  of  the  LLNL  value  (not  shown);  this  is  considered  good  accuracy.  Thus,  the  overall  model  accuracy  is 
not  dependent  on  the  Cp  and  hi  two  thermal  property  assignments  (ours  versus  LLNL),  and  instead  depends 
on  the  accuracy  in  the  values  of  Ri. 
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Using  index  k  to  denote  members  of  the  constituent  radical  set,  the  total  constituent  molar  density 

Nc  =  J2Nk,  (11) 

k 

and  constituent  mole  fractions 

Xk  =  Nk/Nc.  (12) 

A  reference  density  for  nitrogen  N2  is  found  from  the  dry  air  value  at  pref  and  Tref,  giving  Nref=  31.5 
mol/m3.  Using  Nref,  a  dimensionless  N2  molar  density  is  defined  as 

N*=  NNjNref  (13) 


which  acts  as  a  convenient  surrogate  variable  for  p  as  the  N*  value  is  essentially  rate  invariant.  That  is, 
the  partial  pressure  of  N2  is  the  overwhelming  contribution  to  p.  Basically,  the  ratio  of  N*  to  po  is  nearly 
constant.  The  net  constituent  rate  Kc  is  defined  by 


dNc 

dt 


-kcnc. 


A  dimensionless  temperature  is  also  defined 


_  T  —  Tin(N*) 


(14) 


(15) 


where  qualitatively  Tin(N*)  is  the  temperature  at  which  Kc  has  a  value  large  enough  to  signify  combustion 
initiation,  i.e.  positive  dT/dt,  and  Ts  is  such  that  Nc  decreases  from  its  initial  value  by  three  orders  of 
magnitude  (delving  into  higher  order  of  magnitude  decrease  runs  the  risk  of  encountering  round-off  and 
truncation  errors)  at  9  =  0.64,  a  value  chosen  so  that  all  9  values  remain  below  unity  for  all  test  case 
calculations.  A  plot  of  T  versus  6  for  stoichiometric  heptane  combustion  is  shown  in  Fig.  2  where  the  spread 
of  values  for  Ts  (the  slope)  and  (the  ordinate  at  the  origin)  are  shown.  The  variation  of  Tin  is  weak. 
As  N*  increases,  the  change  in  slope  decreases.  Showing  the  appropriateness  of  the  scaling  with  9  rather 
than  using  T,  the  evolution  of  Nc  scaled  by  (po  x  <f>)  for  p0  =  1  bar  and  10  bar  and  <f>  =  0.5,  1.0  and  2.0 
is  illustrated  in  in  Fig.  3.  Following  the  rule  for  finding  Ts  stated  above,  for  <j>  =  0.5  and  2.0  its  value  is  a 
factor  of  0.67  and  1.17,  respectively,  that  of  Ts  at  (f>  =  1.0.  The  normalization  attainable  with  scaling  is  clear. 
Although  it  would  be  tempting  to  curve  fit  dNc/dt  from  Fig.  3  and  adopt  it  as  a  representative  kinetic  rate, 
the  fact  is  that  the  scaling  hides  the  precise  detailed  formulation  that  is  required  for  a  sufficiently  accurate 
kinetics  to  predict  T. 

For  stoichiometric  combustion,  a  plot  of  KC(9,N*)  is  depicted  in  Fig.  4.  Quantitatively,  the  value  of  T)ra 
corresponds  to  a  value  of  Kc  that  is  2%  of  the  minor  first  peak  in  Kc,  with  Tin  =  0(800  K).  The  value  of 
I(c  in  Fig.  4  represents  in  fact  only  the  small  difference  between  the  much  larger  production  and  loss  rates 
which  are  nearly  equal  (not  shown).  This  rate  difference  is  about  1.5%  to  2%  of  either  production  or  loss 
rate  at  the  minor  peak,  reaching  a  maximum  at  9  ~  0.1,  just  after  the  minima  in  Kc ,  then  decreasing  for 
9  >  0.2  to  lows  of  0(0.1%)  of  the  nearly  equal  production  or  loss  rates.  The  maximum  difference  between 
production  and  loss  rates  is  about  5%  for  N*  <  1,  then  decreases  to  about  3%  at  larger  N* .  These  maximum 
differences  tend  to  occur  during  very  large  values  of  dT / dt.  The  small  difference  between  production  and  loss 
is  accurately  captured,  as  the  error  tolerance  in  time  marching  is  10~9,  resulting  in  approximately  a  7  place 
number  accuracy.  The  relatively  small  difference  between  Nc  production  and  loss  suggests  that  Kc  could  be 
treated  as  a  rate  of  a  nearly  quasi-steady  process.  Additionally,  the  region  9  <  0.15  is  where  T  changes  are 
relatively  slow  except  near  9  =  0.15,  but  the  degree  of  quasi-steadiness  based  on  the  production/loss  rate 
differences  is  also  lowest.  The  9  <  0.15  range  is  a  combustion  incubation  region,  with  reaction  times  values 
in  ms  units.  The  9  >  0.2  range  is  the  fast  reaction  region  with  reaction  times  values  in  ps  units.  Curve  fits  to 
Kc  in  the  two  regions  are  plotted  in  Fig.  5.  Shown  on  each  Fig.  5a  and  Fig.  5b  are  5  decades  of  kinetic  rates 
and  the  agreement  on  the  logarithmic  scale  is  very  good.  Each  test  calculation,  i.e.  given  N* ,  is  represented 
in  Fig.  5  by  a  symbol  type;  the  symbol  location  represents  a  selected  sample  of  time  stations  during  the 
calculation.  The  curve  fits  are  based  on  all  time  stations  at  fixed  intervals  of  A 9.  The  first  interval  is  used  for 
9  <  0.06  and  is  small  (~  10~3);  the  second  interval  is  5  times  larger  and  is  used  for  the  remaining  9  values. 
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Mean  heat  capacities  and  enthalpies,  calculated  as  averages  over  the  heavy  (i.e.  decomposed)  species 
using  the  full  rate  (i.e.  LLNL)  model  for  JVj  and  Ri  along  with  Tables  1  through  3  values,  are  defined  as 


Cpm  — 


V 


zC(heavy  species) 


CPtiNi 


H„ 


E 

iC(heavy  species) 


KCNC ' 


(16) 


Plots  of  these  means  are  presented  in  Fig.  6.  Clearly,  the  extrema  in  Hrn  as  function  of  9  occur  at  locations 
that  are  nearly  independent  of  N*;  this  behavior  does  not  occur  when  scaling  by  T  rather  than  9.  In 
particular,  Hm  is  only  weakly  dependent  on  po  (i.e.  IV*),  featuring  an  apparent  near  singularity  as  9  — >  0. 
The  steep  gradients  in  Hm  and  Kc  near  9  =  0  indicate  that,  because  these  Hm  and  Kc  values  are  multiplied 
in  the  energy  equation,  considerable  care  must  be  devoted  to  accurate  modeling  of  the  source  terms  in  this 
incubation  region.  The  ‘noise’  in  Fig.  6a  for  9  >  0.6  is  due  to  numerical  artifacts  as  Nc  — >  0.  The  variation 
of  Cpm  shown  in  Fig.  6b  is  rather  modest,  which  means  that  steep  Hm  gradients  are  due  to  the  heat  of 
combustion  rather  than  to  the  sensible  enthalpy  contribution.  This  observation  focusses  the  models  on  the 
Rk  value  and  evolution,  rather  than  on  the  direct  T  contribution  to  the  Hm  value. 

So  far,  it  is  the  kinetic  rate  of  all  the  constituents,  Kc  that  was  fitted  (Figs.  5a  and  5b).  The  interest 
is  now  in  examining  the  mole  fraction  of  the  constituents,  Xk,  to  determine  their  behavior.  If  the  situation 
were  an  equilibrium,  then  Xk  would  only  depend  on  ( N*,9 ,  (j)).  Since  the  quasi-steadiness  of  Nc  evokes  small 
departures  from  equilibrium,  it  is  natural  to  consider  their  timewise  evolution  through  changes  versus  0; 
thus  Xk  become  functionals  of  9  and  N*  at  fixed  4>.  For  stoichiometric  combustion,  plots  of  8  dominant 
Xk  are  presented  in  Fig.  7.  The  evolution  of  these  constituents  was  obtained  with  the  full  (i.e.  LLNL) 
kinetics  (implicitly  including  thermal  properties)  and  serves  to  indicate  what  are  the  necessary  attributes  of 
the  reduced  model.  The  dependency  on  N*  is  only  pronounced  as  Nc  — >  0,  i.e.  for  large  9.  The  constituent 
radicals  HOO,  00  and  C2H2  reach  peak  Xk  values  of  O(10~2)  in  the  incubation  region,  then  rapidly  decrease. 
The  peak  occurs  at  9  «  0.025  for  the  first  two  radicals,  and  at  9  ~  0.1  for  the  last  radical.  Similar  to  Fig.  6, 
extrema  tend  to  occur  at  fixed  9  values  confirming  again  the  appropriateness  of  the  T  scaling  to  9.  For  further 
modeling  purposes,  it  may  be  noted  that  the  partial  heats  of  formation  for  HOO,  OO  and  O  stated  in  Table 
1  depend  on  na,  the  number  of  carbon  atoms  in  the  radicals  that  are  attached  to  the  particular  constituent 
radical.  Because  the  variation  of  Sh°  with  species  identity  for  those  particular  constituents,  average  values 
are  needed  to  calculate  the  T  variation  from  the  energy  equation.  The  average  values  calculated  using  the 
LLNL  model  are  -113  kJ/rnol  for  HOO  and  O,  and  -42kJ/mol  for  OO. 

Aside  from  air  and  combustion  products,  the  light  species  CO,  H2,  CH4,  H2O2,  CH20,  C2H4  and  C2H2 
have  the  largest  molar  density  value  of  the  light  species  in  the  base.  The  CO  molar  density  dominates  that 
of  each  of  the  other  species  in  the  list. 


III.  Summary  and  Conclusions 

The  species  entering  oxidation  kinetics  were  categorized  into  either  heavy  or  light  species.  For  heavy 
species  (i.e.  those  with  carbon  number  >  3),  a  constituent  radical  decomposition  has  been  proposed.  The 
constituent  radicals  and  the  light  species  form  a  base  of  small  dimension  used  to  follow  the  kinetics.  It  has 
been  shown  that  modeled  simplified  expressions  for  constituent  enthalpies  lead  to  an  adequately  modeled 
heat  release  for  heptane  combustion.  A  simplified,  low-dimensional  chemistry  model  thus  depends  on  an 
adequate  representation  of  a  reduced  rate  set,  Rk,  of  the  constituents.  The  net  constituent  rate,  Kc,  is 
nearly  quasi-steady  and  may  be  split  into  an  incubation  region  of  modest  temperature  rise  followed  by  a  fast 
reaction  region  of  high  temperature.  A  degree  of  normalization  of  molar  density  behavior  is  achieved  using 
scaled  pressure  and  temperature  variables,  N*  and  9,  respectively. 

The  incubation  (small  9)  region  is  characterized  by  roughly  ms  reaction  times.  These  times  are  similar 
to  diffusion  time  scales,  and  thus  these  processes  are  expected  to  significantly  couple  with  the  chemistry 
during  incubation.  When  reducing  the  species  to  constituents,  heavy  species  diffusion  must  be  represented 
in  this  new  framework.  The  question  is  then:  What  are  appropriate  diffusion  coefficients?  Since  a  particular 
constituent  spans  the  entire  heavy  species  set,  a  heavy  species  mean  diffusion  coefficient  is  appropriate. 
This  means  that  heavy  species  differential  diffusion  cannot  be  considered  in  the  incubation  region.  This  is 
a  problem  encountered  in  any  reduction  scheme  that  focusses  on  modeling  both  the  fully  developed  flame 
and  the  incubation  region  with  the  same  reduced  rate  set.  On  the  other  hand,  hypothetically,  if  a  species  is 
crucially  affecting  the  incubation  kinetics  through  differential  diffusion,  it  can  be  selected  to  be  part  of  the 
reduced  set.  Generally,  using  a  reduced  model  in  the  incubation  region,  means  that  accuracy  in  this  initial 
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evolution  path  is  expected  to  be  somewhat  limited.  This  has  implications  for  predicting  ignition  delay  times. 

In  contrast  to  the  incubation  region,  in  the  fast  reaction  zone  the  coupling  between  chemistry  and  flow 
processes  is  weak  and  combustion  is  primarily  determined  by  the  mixing  rate.  The  behavior  in  the  test  cases 
using  the  LLNL  species  and  rate  data  is  that  of  very  rapid  temperature  rise  beginning  at  9  «  0.1,  i.e.  at 
the  transition  of  the  incubation  region  to  that  of  fast  rate  reactions.  The  temperature  profiles  in  the  fast 
rate  region  tend  to  be  independent  of  the  details  of  behavior  during  incubation,  that  is,  an  approach  to  a 
high-temperature  ‘attractor  state’  seems  to  prevail.  Fortunately,  this  implies  that  any  inaccuracies  developed 
during  incubation  will  be  eventually  eliminated.  However,  as  mentioned  above,  special  accuracy  in  modeling 
of  the  incubation  region  is  needed. 

Model  closure  will  require  the  determination  of  effective  mean  source  strengths  for  the  light  species  and 
light  radicals  resulting  from  the  heavy  species.  The  final  model  will  thus  focus  on  reactions  of  the  light 
species;  the  necessary  rates  are  expected  to  be  well  determined  insofar  as  kinetic  interactions  among  light 
species  prevail. 

The  constraint  of  atomic  conservation  enables  a  reduction  in  the  number  of  rate  progress  variables  by  the 
number  of  atoms  in  the  system.  Formally,  this  could  be  applied  to  air  and  combustion  products.  However, 
from  the  viewpoint  of  numerical  accuracy,  these  constraints  are  best  applied  to  species  (other  than  N2),  or 
species  combinations,  of  relatively  small  mole  fraction  in  order  to  minimize  effects  of  round-off  errors. 
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Radical 

Sh° 

5hc 

5hf 

ah 

bh 

CH3  (methyl) 

-42.0 

714 

188 

3.137 

3.433 

CH2  (methylene) 

-20.8 

614.3 

411 

2.784 

2.812 

C'H  (methylidyne) 

£*  -7 

^  507 

^  603 

2.068 

3.038 

C2H3  (vinyl) 

62.5 

1212 

237 

4.432 

4.468 

C2H2  (CH=CH)® 

73  ±  2 

1102  ±  2 

not  relevant 

3.82 

3.72 

C2  (C=C,  dicarbon) 

230 

1017 

608 

2.919 

1.796 

HC2  (HCeeC,  ethynyl) 

227 

1135 

339** 

3.937 

2.2515 

CO  (keto) 

-133  ±  2 

260.5 

23 

2.892 

1.531 

HCO  (formyl) 

-124  ±  3 

390.3 

167 

3.5 

2.1 

HO  (hydroxyl) 

-166  ±  6*** 

-45 

204 

2.084 

1.0945 

HOO  (hydroperoxy)*® 

-205+116/na 

-84+116/na 

215-116/n0 

4.194 

1.908 

00  (peroxy)* 

-16-13na 

-16-13na 

not  relevant 

3.65 

2.135 

O  (ethers)*# 

-74-13na 

-74-13na 

323+13na 

-1.46 

2.44 

®trans:  -2,  cis:  +2;  for  a  molecule  with  one  or  more  vinyl  also  present  (i.e.,  polyene), 

SCp/Ru  is  increased  by  12%; 

** Alternate  value  (op.  cit.  NIST)  is  250  from  older  data; 

***value  =  -180  if  combined  with  CH  radical; 

*na  is  the  number  of  carbon  atoms  in  attached  radicals; 

®effective  na  =  1.16  with  ring  or  branched  attached  radical; 

#each  of  the  two  attachment  points  has  a  contribution  to  na,  A na,  which  has  a  maximum  value  of  2. 

The  effective  na  =  3  for  a  ring  with  na  >  4. 

If  at  least  one  vinyl  is  also  present,  then  add  the  increment  Ana  =  +1. 

For  a  ring  with  na  =  2  or  3,  a  correction  of  90  must  be  added  to  the  Sh°  value  listed  in  the  table. 
5CP/RU  fit  is  for  ring  ethers;  for  a  chain,  6CP/RU  =  2.5. 


Note:  if  the  keto  (or  formyl)  radical  is  attached  to  oxygen  (or  hydroxyl),  then  a  correction  of  -90  must  be 

added  to  the  Sh°  value  listed  in  the  table. 

Table  1.  Thermodynamic  properties  of  constituent  radicals.  5h° ,  5hc,  ShJ  heats  of  formation,  combustion,  and 
component-to-free  transition,  respectively,  in  kJ/mol;  constants  ar  and  bh  for  partial  molar  heat  capacity  in 
the  form  5CP/RU  =  ah  +  bhln(T/Tref)\  Tref  =  25  C  =  298.15  K. 
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Radical 

Sh*  (kJ/mol) 

ah 

bh 

alkyl 

182  (174)* 

0.6 

-1.5 

alkenyl 

235  (128)* 

0.0 

0.0 

alkynyl 

240  (133)* 

1.2 

-1.5 

Q  +  OOH** 

239 

0.6 

-2.8 

alkyl  +00** 

78 

-0.7 

-1.4 

alkyl  +  O 

147# 

-1.3 

0.5 

alkyl  +  CO 

160# 

0.2 

-0.6 

Q  +  00  +  OOH 

140  (70)* 

1.3 

-3.8 

Ring  +  O 

226 

3 

-4 

dienes  (e.g.  C4H5) 

196  (90)* 

0.0 

0.0 

Stable  species  Shc  (kJ/mol)  ah  bh 


Ring  ether  90  (80  with  2  CH)  2.2  -2.5 

nc[n]ket  not  relevant  1.3  -1.8 

*  value  with  CH  constituent  or  single  bond  C; 

**for  na  ^  3;  if  OOH  is  on  a  branch,  use  alkyl  values; 

#with  alkenyl,  subtract  20. 

Table  2.  Rules  for  enthalpy  corrections  for  free  radicals  and  enthalpy  corrections  for  a  few  stable  heavy 
molecules.  SCP/RU  =  ah  +  bhln(T/Tref)  for  heavy  free  radicals  and  also  corrections  for  ring  ethers  and  ketohy- 
droperoxy’s  (only  nc[n]ket,  where  n  is  the  carbon  count  in  the  species). 
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Molecule /Radical 

h° 

hc 

ah 

bh 

h2 

0.0 

241.5 

3.282 

0.400 

o2 

0.0 

0.0 

3.476 

0.5663 

n2 

0.0 

0.0 

3.388 

0.469 

c 

717 

1111 

2.50 

0.0 

h2o 

-241.5 

0.0 

3.688 

1.217 

co2 

-393.5 

0.0 

4.690 

1.390 

N 

473 

473 

2.50 

0.0 

H 

218.0 

339 

2.50 

0.0 

HO 

38 

159 

3.385 

0.3637 

H00 

10.5 

131 

4.150 

1.307 

0 

249.2 

249.2 

2.536 

0.0 

CO 

-110.5 

283 

3.426 

0.4749 

HCO 

43.1 

558 

4.154 

1.2875 

ch4 

-74.6 

802 

3.797 

4.305 

ch3 

146 

902 

4.440 

2.249 

ch2 

390 

1025 

3.973 

1.3015 

CH 

596 

1110 

3.220 

0.7136 

c2h3 

300 

1449 

5.1 

3.5 

hc2 

566 

1474 

4.434 

1.404 

C2 

838* 

1625 

4.58 

0.0 

NO 

90.3 

90.3 

3.533 

0.4508 

no2 

33.2 

33.2 

4.691 

1.151 

h2o2 

-136 

106 

5.269 

1.880 

HCOH 

-109 

526.5 

4.27 

2.546 

c2h2 

228 

1257 

5.368 

2.294 

c2h4 

52.5 

1323 

5.383 

4.676 

*  Alternate  value  of  832  from  the  CRC  tables. 


Table  3.  Thermodynamic  properties  of  molecules  and  free  radicals.  h°  and  hc  (heats  of  formation  and  combus¬ 
tion,  respectively),  in  kJ/mol;  constants  ah  and  bh  for  molar  heat  capacity  in  the  form  Cp/Ru  =  ah  +bhln(T /Tref)-, 
Tref  =  25C  =  298.15  K. 
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cc 

a 

o 


(c) 

Figure  1.  Cp/Ru  versus  log(T)  where  T  is  in  degrees  K.  Data  values  are  in  symbols;  lines  portray  the  fit 
constants,  (a)  alkanes  (constants  from  Table  1);  (b)  alkenes  (constants  from  Table  1);  (c)  ring  ethers  (CH2)nO 
(constants  from  Table  1);  and  (d)  some  free  radicals  (constants  from  Table  3). 
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Figure  2.  Temperature  (in  K)  versus  6  for  various  values  of  N* .  po  is  in  bar. 


Figure  3.  Evolution  of  the  normalized  Nc  versus  6  for  po  =  1  bar  and  10  bar  for  <f>  =  0.5,  1.0  and  2.0. 
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Figure  4.  Kc(0,  N*)  versus  6  for  stoichiometric  combustion  of  heptane  from  test  runs,  po  is  in  bar. 


Figure  5.  Kc  fits  of  the  test  runs  in  Fig.  4  for  (a)  the  incubation  region,  and  (b)  the  fast  reaction  region. 
Curves  denote  fits  and  the  legend  is  the  same  as  in  Fig.  4;  symbols  denote  selected  full  rate  model  values  (the 
fits  were  based  on  all,  i.e.  not  only  the  shown  selected,  values). 


Figure  6.  Mean  values  over  the  heavy  species  set  as  functions  of  6  and  N* .  (a)  Hm/(RuTref)  and  (b)  CPm/Ru- 
po  is  in  bar. 
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Figure  7.  Constituent  mole  fractions  X &  as  functions  of  6  and  N*  for  0=1.  The  constituent  is  (a)  CH2;  (b) 
CH3;  (c)  C2H3;  (d)  CH;  (e)  CO;  (f)  HCO;  (g)  HO;  and  (h)  O.  p0  is  in  bar. 
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A  simplified  model  is  proposed  for  the  kinetics  of  alkane  oxidation  in  air,  based  on  a 
decomposition  of  heavy  (carbon  number  greater  or  equal  to  3)  hydrocarbons  into  a  13 
constituent  radical  base.  The  behavior  of  this  base  is  examined  in  test  computations  for 
n-heptane  utilizing  Chemkin  II  with  LLNL  data  inputs,  placing  emphasis  on  modeling  to 
predict  the  heat  release  and  temperature  evolution.  A  normalized  temperature  was  con¬ 
structed  which  when  used  to  plot  the  total  constituent  molar  density  divided  by  the  product 
of  the  equivalence  ratio  and  a  nondimensional  pressure,  reveals  a  self-similar  behavior  of  the 
plotted  variable  over  a  wide  range  of  initial  pressures  and  equivalence  ratios.  Examination 
of  the  LLNL  kinetics  shows  that  the  total  constituent  molar  density  rate  follows  a  quasi¬ 
steady  behavior.  This  reaction  rate  was  curve  fitted  along  with  the  corresponding  enthalpy 
production.  The  fits  are  shown  against  the  normalized  temperature  for  various  equiva¬ 
lence  ratios  and  initial  nondimensional  pressures  and  comparisons  with  the  LLNL  kinetics 
are  very  favorable.  The  model  reduces  the  LLNL  n-heptane  mechanism  from  160  species 
(progress  variables)  and  1540  reactions  to  12  progress  variables,  16  quasi-steady  rates  (as¬ 
sociated  with  heavy  species),  162  conventional  reaction  rates  (light  species)  and  11  other 
functional  forms  (i.e.  fits  for  the  mean  heavy-species  heat  capacity  at  constant  pressure, 
the  enthalpy  release  rate  of  the  heavy  species,  and  the  molar  fraction  of  quasi-steady  light 
species).  The  proposed  kinetic  mechanism  is  valid  over  a  pressure  range  from  atmospheric 
to  60  bar,  temperatures  from  600  K  to  2500  K  and  equivalence  ratios  from  0.125  to  8.  This 
range  encompasses  diesel,  HCCI  and  gas  turbine  engines,  including  cold  ignition;  and  NOj, 

CO  and  soot  pollutant  formation  in  the  lean  and  rich  regimes,  respectively. 

I.  Introduction 

Computational  models  for  turbulent  reactive  flows  present  challenges  that  have  yet  to  be  fully  met. 
First,  combustion  of  typical  hydrocarbons  proceeds  through  the  formation  of  hundreds  of  species  interacting 
through  thousands  of  reactions.  Second,  a  complete  numerical  resolution  of  high  Reynolds  number  flow 
requires  a  very  large  number  of  grid  points  that  is  beyond  current  typical  supercomputer  capabilities.  Thus, 
combining  detailed  reaction  rate  models  with  the  fundamental  flow  equations  yields  a  numerical  problem 
that  is  computationally  prohibitive.  A  drastic  (at  least  order  of  magnitude)  reduction  in  the  dimension  of 
reaction  kinetics  is  required  for  a  feasible  computational  model.  In  spirit,  this  reduction  is  similar  to  that 
required  in  Large  Eddy  Simulations  (where  energetically  significant  flow  scales  are  computed  and  the  others 
are  modeled)  with  respect  to  Direct  Numerical  Simulations  (where  all  flow  scales  are  computed).  Since  each 
species  has  a  characteristic  time  scale,  in  the  kinetic  reduction  it  is  also  desirable  to  compute  only  those 
species  having  an  essential  characteristic  time  scale  (to  be  defined)  and  model  the  kinetics  of  the  remaining 
species. 

We  categorize  all  species  involved  in  hydrocarbon  oxidation  into  light  species  and  heavy  species;  the  heavy 
species  are  those  having  a  carbon  number  n  ^  3.  Based  on  this  definition,  a  method  for  simplifying  alkane 
combustion  is  proposed  and  tested  for  n-heptane.  The  method  focusses  on  the  total  constituent  radical  molar 
density  (see  Section  II  definition  and  constituent  radical  base  formation)  of  the  heavy  species  rather  than 
on  the  much  larger  number  of  these  species.  That  is,  rather  than  following  all  species  through  their  reaction 
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coordinates,  we  follow  a  reduced  set  of  reaction  coordinates;  this  reduced  set  is  called  a  base.  Our  primary 
objective  is  the  accurate  determination  of  the  energetics,  i.e.  heat  release  and  temperature  evolution;  at  this 
initial  stage  of  model  development,  neither  NOa,  nor  soot  kinetics  are  addressed.  In  Section  III,  the  behavior 
of  the  base  set  is  examined  in  the  context  of  a  detailed  reaction  model  for  heptane  combustion  as  given 
by  LLNL1  and  the  proposed  model  is  tested  against  the  LLNL  kinetic  scheme.  Section  IV  is  devoted  to  a 
discussion  on  the  completion  of  the  model.  Finally,  in  Section  V  a  summary  of  the  accomplishments  and  a 
perspective  on  the  future  work  are  offered. 


II.  Species  representation:  a  reduced  base  set 


A.  Heavy  species  decomposition  into  constituent  radicals 

Consider  the  count  of  atoms  (H,C,0,N)  in  species  i.  We  define  Wji  as  the  number  of  species  j  in  species  i 
that  produces  the  correct  atom  count,  where  the  index  j  denotes  the  species  set  of  air  (02,N2)  and  final 
combustion  products  (H20,C02).  (This  means  that  mathematically,  we  consider  species  i  as  being  composed 
of  air  and  final  combustion  product  species.)  The  molar  densities  of  species  in  this  set  are 


Nj=N9  (!) 

where  N®  is  the  effective  total  molar  density  of  species  j.  The  same  exercise  of  atom  counting  and  species 
representation  by  the  j  set  can  be  performed  for  all  species  in  the  mixture  to  obtain  N®.  That  is,  the  atomic 
elements  molar  densities  are  represented  by  species  j  molar  densities.  Conservation  of  atoms  in  the  reactions 
makes  unnecessary  utilizing  rate  expressions  for  reaction  coordinates  of  air  species  and  complete-combustion 
products,  as  a  one-to-one  map  exists  between  these  and  the  atom  set.  The  heavy  species  constituent  radicals 
(CH2,  CH3,  CH,  C2H3,  C2H2,  C2,  HC2,  CO  (keto),  HCO,  HO,  H02,  00,  O)  form  a  set  of  13  entities  from 
which  any  heavy  species  or  radical  may  be  constructed.  The  total  constituent  molar  density  is 

13 

Nc  =  ^Nk  (2) 

k= 1 


where  the  molar  density  of  constituent  k  is  Nk-  The  decomposition  into  constituents  is  similar  to  group 
additivity,2, 3  except  that  it  considers  only  first  order  (compositional)  effects  and  only  interactions  with 
adjacent  groups.  The  reaction  rate  associated  with  Nc  is  defined  by 

K'  = - IS—  ,3) 

We  define  a  reference  density  for  nitrogen  N2  from  the  dry  air  value  at  pressure  pref  =  1  bar  and  Tref= 
298.15  K,  giving  Nref=  31.5  mol/m3.  Using  Nref,  a  dimensionless  N2  molar  density  is  defined  as 

N*  =  NN2/Nref  (4) 


which  acts  as  a  convenient  surrogate  variable  for  the  pressure  p  as  the  N*  value  is  essentially  rate  invariant. 
That  is,  the  partial  pressure  of  N2  is  the  overwhelming  contribution  to  p.  Basically,  the  ratio  of  N*  to  po  is 
nearly  constant.  A  normalized  temperature  variable  6  is  defined  by 


T-Ts 

Tr(<P,  N* 


Tr  =  2065(1V*)U  «#),  w(4>)  =  0 


1.5  +  1.310 
1  +  0.710+1.1  02 


(5) 


where  Ts  is  a  modeling  parameter  representing  the  lowest  temperature  at  which  Kc  >  0.  The  9  definition  is 
such  that  Nc  decreases  during  stoichiometric  combustion  from  its  initial  value  by  three  orders  of  magnitude 
(delving  into  higher  order  of  magnitude  decrease  runs  the  risk  of  encountering  round-off  and  truncation 
errors)  for  6  >  0.6,  a  value  chosen  so  that  all  9  values  remain  below  unity  for  all  test  case  calculations. 

The  contribution  of  the  heavy  group  to  the  rate  of  change  of  a  light  species  molar  density  A*  with  mole 
fraction  X{  (relative  to  the  light  group)  is  expressed  in  terms  of  quasi-steady  gain  and  loss  rates  K Gj  and 


KLi  as 


dNi 


dt 


heavies 


Nc(KGi  -  XiKLi). 


(6) 
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Aside  from  the  molar  reaction  rates,  determining  the  temperature  evolution  in  a  reactive  system  requires 
knowledge  of  the  species  molar  enthalpies  and  heat  capacities.  For  species  i,  the  molar  enthalpy  may  be 
expressed  as  ht  =  h®  +  hi(T)  where  h(-  is  the  heat  of  formation,  ht  is  the  sensible  enthalpy  and  T  is  the 
temperature.  The  heat  of  formation  is  taken  at  the  above  reference  conditions,  giving  hi  =  flj  ^  Cp^dT 
where  CPti  is  the  molar  constant-pressure  heat  capacity.  A  heat  of  combustion  for  species  i  is  given  by 

K  =  hi  ~YlW^h°3  (7) 

3 

under  the  assumption  that  water  remains  in  vapor  state.  The  energy  equation  is 


NrQ 


p,h 


^  TV  T  \  dT 

^  ]  Cp^Ni  |  -7- 
i(zlights 


dt 


hiRi  +  Nc(RuTref)Kh 

i(zlights 


where 

'Ri  =  {dNi/  dt)reac, 

Ru  is  the  universal  gas  constant, 

Cp,h  =  (  ^  CPtlNt)/Nc, 

l(zheavies 

and  the  rate  of  enthalpy  change  due  to  heavy  species  l  chemical  kinetic  changes  is  defined  as 


Kh  = 


hiRi 


1 


KlEheavies 


R„TrefN, 


(8) 

(9) 

(10) 


(11) 


B.  Light  species  as  part  of  the  base 

Light  species  are  not  subject  to  meaningful  decomposition  according  to  the  above  procedure.  Examination 
of  runs  made  using  the  LLNL  database  indicate  that  the  light  species  should  be  categorized  in  two  subsets. 
The  species  in  the  first  subset  require  rate  equations  and  are  part  of  the  base  set.  The  species  in  the  second 
subset  have  a  quasi-steady  behavior  and  require  fits  of  their  mole  fraction  value  as  a  function  of  the  state 
variables  (T,N*,<j>)  and  modeling  parameters  (here,  only  Ts).  This  first  subset  consists  of:  H20,  C02,  02, 
H,  CO,  H2,  CH4,  H202,  C2H2i  C2H4,  CH20.  The  second  subset  are  the  radicals:  O,  CH,  CH2,  CH3,  HO, 
HCO,  H02,  HC2,  C2H3.  Thus  this  entire  set  consists  of  20  light  species.  For  their  energetics,  values  of  hc 
are  fitted  from  the  literature  and  Cp  s  are  obtained  as 

Cp/Ru  =  ah  +  bh\n(T/Tref).  (12) 

Other  properties  of  the  light  molecules  and  light  free- form  of  the  constituent  radicals  are  given  in  Table  l.4  7 

C.  Final  base  set 

The  reaction  base  set  is  composed  of  a  global  (i.e.  total)  constituent  radical  and  11  light  molecules  or  free 
radicals  of  the  first  subset  (Table  1).  Thus,  there  are  only  12  progress  variables.  To  compute  Nc  one  must 
model  Kc ;  to  compute  the  11  light  species,  one  must  model  the  KGi,  the  KLi,  the  Xi  of  the  second  subset 
and  the  conventional  light  species  reaction  rates  of  the  first  subset. 

III.  Results 

A.  Identification  of  a  self-similar  behavior 

Showing  the  appropriateness  of  the  scaling  with  9  rather  than  using  T,  the  evolution  of  Nc  scaled  by  (<j>  x  N*) 
for  fixed  po  at  several  <f>  values  and  at  fixed  </>  for  several  po  values  is  illustrated  in  Fig.  1  using  the  LLNL 
kinetics.  The  normalization  attainable  with  9  scaling  is  clear  and  so  is  the  fact  that  Nc/(<j>  x  N*)  is  a 
similarity  variable.  At  fixed  po  =  20  bar,  Fig.  la,  all  Nc/(<j>  x  N*)  curves  coincide  nearly  on  the  same 
curve,  with  the  somewhat  dissenting  curve  for  <f>  =  4  showing  that  for  very  rich  mixtures  there  are  indeed 
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constituents  left  after  the  oxidizer  was  depleted.  However,  even  for  as  rich  a  mixture  as  (f>  =  2,  similarity 
holds,  making  this  model  valid  in  realms  beyond  those  in  advanced  reduced  schemes8  where  the  upper  limit 
of  the  scheme  validity  is  </>  =  1.5,  at  most9  (j>  =  2.0  or  exceptionally10  (f>  =  3.  The  behavior  in  Fig.  la  is 
typical  of  all  fixed  po  plots.  In  Fig.  lb,  corresponding  plots  at  stoichiometric  conditions  are  shown,  which  are 
generic  of  all  fixed  </>  plots.  Again,  the  similarity  achieved  with  the  two  normalized  variables  is  remarkable, 
and  only  at  a  combination  of  high  po  values  and  small  Ts  values  (typical  of  the  lowest  T  limit  at  which  Kc 
rates  are  finite),  are  small  to  moderate  departures  from  the  excellent  normalization  perceptible. 

We  note  that  by  9  ~  [0.55, 0.6],  Nc  0,  meaning 
that  modeling  errors  past  that  9  range  are  somewhat 
irrelevant  to  the  model’s  accuracy  since  the  reaction 
is  practically  finished.  Although  it  would  be  tempt¬ 
ing  to  curve  lit  dNJdt  from  Fig.  1  and  adopt  it  as  a 
representative  kinetic  rate,  the  fact  is  that  the  scal¬ 
ing  hides  the  precise  detailed  formulation  that  is  re¬ 
quired  for  a  sufficiently  accurate  kinetics  to  predict 
T.  Thus,  fits  for  Kc  are  instead  needed  to  achieve 
appropriate  accuracy,  as  discussed  below. 

B.  Global  constituent  molar  density  model 

Plots  of  Kc  versus  9  show  that  there  are  two  differ¬ 
ent  regions  of  the  curve.  For  small  9  values,  there  is 
a  combustion  incubation  region,  with  reaction  times 
values  in  ms  units.  As  9  increases,  one  encounters 
the  fast  reaction  region  with  reaction  times  values 
in  p s  units.  Numerous  tests  made  with  the  LLNL 
kinetics  showed  that  Kc  represents  only  a  small  dif¬ 
ference  between  production  and  loss  rates  which  are 
nearly  equal  (not  shown)  and  are  much  larger  than 
Kc.  The  Kc  value  is  typically  about  1.5%  to  5%  of 
either  production  or  loss  rate.  The  maximum  differ¬ 
ence  between  production  and  loss  rates  is  about  5% 
for  N*  <  1,  then  decreases  to  about  3%  at  larger 
N*.  These  maximum  differences  tend  to  occur  dur¬ 
ing  very  large  values  of  dT/dt.  This  fact  suggests 
that  Kc  could  be  treated  as  a  rate  of  a  nearly  quasi¬ 
steady  process.  If  fits  of  Kc  can  be  obtained  over  a 
significant  region  of  parametric  behavior,  this  model 
would  be  compact  and  represent  a  considerable  re¬ 
duction  in  the  problem  dimensionality.  Noteworthy, 
had  the  approach  been  to  find  Kc  from  slopes  of  Nc 
in  Fig.  1,  it  is  clear  that  the  important  incubation 
region  would  have  been  missed. 

Fitting  Kc  as  a  function  of  T  and  parameters 
N*,  (f)  and  Ts  presented  a  major  challenge  because 
of  the  non-monotonic  Kc  behavior  and  also  because 
Kc  spans  so  many  decades.  Fits  were  obtained  by 
first  inspecting  the  LLNL  kinetics  to  understand  the 
rate  constant,  K ,  variation.  This  variation  showed 
an  incubation  region  (9  <  0.2)  featuring  first  a  max¬ 
imum,  Kmx,  followed  by  a  minimum,  A'm„,  and  a  high-T  region  ( 9  >  0.2)  where  Kc  monotonically  increases. 
Thus,  the  functional  fits  are  performed  in  three  separate  regions,  with  connection  constraints  between  con¬ 
secutive  regions  (the  ultimate  fits  are  piecewise  continuous) .  The  non-monotonic  incubation-region  behavior 
is  captured  by  using  a  cubic  functional  mapping  from  K(T)  to  y(T )  through  the  form 


Mo/Ra 

h° 

hc 

ah 

bh 

h2 

0.0 

241.5 

3.282 

0.400 

o2 

0.0 

0.0 

3.476 

0.5663 

n2 

0.0 

0.0 

3.388 

0.469 

c 

717 

1111 

2.50 

0.0 

h2o 

-241.5 

0.0 

3.688 

1.217 

co2 

-393.5 

0.0 

4.690 

1.390 

N 

473 

473 

2.50 

0.0 

H 

218.0 

339 

2.50 

0.0 

HO 

38 

159 

3.385 

0.3637 

HOO 

10.5 

131 

4.150 

1.307 

O 

249.2 

249.2 

2.536 

0.0 

CO 

-110.5 

283 

3.426 

0.4749 

HCO 

43.1 

558 

4.154 

1.2875 

ch4 

-74.6 

802 

3.797 

4.305 

ch3 

146 

902 

4.440 

2.249 

ch2 

390 

1025 

3.973 

1.3015 

CH 

596 

1110 

3.220 

0.7136 

c2h3 

300 

1449 

5.1 

3.5 

hc2 

566 

1474 

4.434 

1.404 

C2 

838* 

1625 

4.58 

0.0 

NO 

90.3 

90.3 

3.533 

0.4508 

N02 

33.2 

33.2 

4.691 

1.151 

h2o2 

-136 

106 

5.269 

1.880 

HCOH 

-109 

526.5 

4.27 

2.546 

C2H2 

228 

1257 

5.368 

2.294 

c2h4 

52.5 

1323 

5.383 

4.676 

*  Alternate  value  of  832  from  the  CRC  tables. 

Table  1.  Thermodynamic  properties  of  molecules  and 
free  radicals.  h°  and  hc  (heats  of  formation  and  com¬ 
bustion,  respectively),  in  kJ/mol;  constants  or  and  bh 
for  molar  heat  capacity  in  the  form  Cp/Ru  =  or  + 
bhln(T/Tref );  Tref  =  298.15  K.  "Mo"  denotes  "mole¬ 
cule".  "Ra"  denotes  "radical". 


K (T)  —  —  ( Kmx  +  Kmn)  +  —(Kmx  ~  Kmn)(y>>  —  3y) 


(13) 
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where  y(T)  =  —  1  at  T  =  Tmx  which  is  the  location  of  Kmx  and  y(T)  =  1  at  T  =  Tmn  which  is  the  location 
of  Kmn.  In  fact,  the  range  of  y  is  from  slightly  smaller  than  -2  at  T  =  Ts  (i.e.  Kc  >  0)  to  y  <  1.8.  The 
values  of  Kmx,Kmn,Tmx  and  Tmn  are  fitted  as  polynomials  in  parameters  ln(7V*),  ln(</>)  and  ln(Ts/Tre/). 
Also  fitted  is  To  which  is  the  T  value  at  which  y  =  0.  Two  different  changes  of  variables  are  made  from  y(T) 
to  y(z)  depending  on  how  T  compares  to  T0.  For  T  <  T0,  a  variable  z  =  (T0  —  T)/(T0  —  Tmx)  is  defined, 
which  means  that  y  =  —  1  at  0  =  1.  To  achieve  the  mapping,  we  first  generate  (y,  z )  pairs  which  are  functions 
of  parameters  ln(iV*),  ln(</>)  and  ln(Ts/Tre/)-  To  match  these  pairs,  a  ratio  of  appropriate  functions  (e.g. 
polynomials  a  +  bz  +  cz2  +  ...;  power  functions  of  type  az;  and  combinations  of  the  two  functional  forms) 
is  chosen  to  produce  the  y(z)  mapping.  For  T  >  To,  i.e.  y  ^  0,  a  similar  procedure  is  used  for  y{z)  where 
now  z  =  (T  —  T0)/(Tmn  —  T0).  Beyond  incubation,  for  the  high-T  region,  ln(A')  is  fitted  as  a  fifth  order 
polynomial  in  9P  where  p  =  1  for  <j>  <  0.5  and  p  =  (0.86  +  0.28^)/(0.74  +  0.52 </>)  for  </>  >  0.5.  Scrutiny  of  the 
LLNL  kinetics  showed  that  the  largest  value  of  Ts  is  825  K.  Then,  for  Ts  =  825  K,  the  polynomial  coefficients 
were  fitted,  separately  for  the  lean/stoichionretric  (</>  ^  1)  and  rich  [<f>  >  1)  regions,  as  polynomials  in  In (N*) 
and  ln(0).  To  obtain  fits  for  Ts  <  825  K,  a  scaling  procedure  was  used  from  the  Ts  =  825  K  fits  to  match 
particular  K  values  at  a  specific  point  in  the  incubation  region. 

All  fits  further  presented  in  Figs.  2-6  have  been  achieved  with  the  fixed  set  of  mapping  functions  de¬ 
scribed  above.  That  is,  no  local  corrections  of  any  sort  were  used,  nor  exceptions  from  those  functions  were 
considered.  Local  corrections  in  a  four-dimensional  space  are  not  computationally  practical. 

Curve  fits  to  Kc  computed  using  the  LLNL  kinetics  are  plotted  in  Figs.  2-4.  In  Fig.  2a  and  2b  (fixed 
Po,  various  <j>  and  Ts),  and  Fig.  3a  and  3b  (fixed  </>  in  the  lean  regime,  various  po  and  Ts)  we  differentiate 
between  the  incubation  and  fast-reaction  regions  discussed  above  and  show  on  each  figure  up  to  5  decades 
of  kinetic  rates.  On  each  figure,  symbols  represent  the  LLNL-based  computation  and  the  symbol  location 
represents  a  selected  sample  of  time  stations  during  the  calculation.  The  curve  fits  are  shown  as  lines  and 
are  based  on  time  stations  at  fixed  intervals  of  A 9.  The  agreement  on  the  logarithmic  scale  is  excellent  to 
very  good  and  deterioration  of  the  fits  at  9  >  0.55  correspond  to  very  small  values  of  Nc,  as  noted  above. 
Particularly  important  is  the  excellent  agreement  between  fits  and  LLNL  kinetics  in  the  region  of  9  <  0.2 
representing  54%  of  the  Nc  decay  (see  Fig.l).  For  9  >  0.2  the  agreement  is  still  excellent  as  long  as  Ts  is 
not  too  low;  good  agreement  still  occurs  in  the  significant  Nc  regions  (i.e.  9  <  [0.55,0.6])  and  over  a  much 
wider  (j>  range  than  usually  shown  in  the  literature.  Figure  4a  focusses  on  the  rich  regime,  showing  the  entire 
range  of  relevant  9  (see  Fig.  la)  at  various  po  and  Ts.  The  agreement  of  fits  with  the  LLNL  data  is  excellent, 
bodying  well  for  future  studies  addressing  soot  and  CO  formation.  Noteworthy,  the  fits  encompass  a  wide 
Po  range  including  the  desired  high-pressure  situations  for  engine  combustion,  and  a  much  wider  (f>  range 
than  typically  used  in  kinetic  reduction  schemes8  where  generally  </>  €  [0.5, 1.5],  sometimes9  <f>  €  [0.5, 2.0] 
and  exceptionally10  (f>  £  [0.7,  3.0].  To  show  the  influence  of  Ts,  in  Fig.  4b  is  depicted  the  case  for  which  the 
widest  range  of  Ts  was  available,  corresponding  to  po  =  10  bar  and  (f>  =  1.  Similar  to  the  other  figures,  the 
agreement  with  the  LLNL  kinetics  is  excellent  to  very  good  up  to  9  <  0.2  and  somewhat  deteriorates  past 
that  station. 

This  deterioration  is  a  manifestation  of  the  need  for  matching  two  different  functional  forms:  one  in  the 
incubation  region  (y  >  0),  and  the  other  in  the  high-T  region.  The  matching  point  is  at  ymax;  as  ymax 
increases,  the  high-T  region  fit  improves  while  the  incubation  fit  worsens,  and  vice  versa  for  decreasing  ymax. 
Fits  can  also  be  improved  if  one  is  less  ambitious  in  covering  the  large  </>  range  treated  here.  The  functional 
behavior  of  (f>  outliers  is  somewhat  different  from  that  for  moderate  <p .  For  </>  <  1  and  small,  the  high-T  rates 
are  quite  modest  and  below  Kmx.  For  </>  >  1,  as  (j>  increases,  the  ratio  Kmx/Kmn  decreases  and  high-T  rates 
exceed  Kmx  except  for  the  largest  <j>  where  the  reaction  terminates  before  the  Kmx  value  is  reached  again  in 
the  high-T  region.  This  behavior  is  independent  of  po-  This  is  to  say  that  if  the  outliers  are  eliminated,  a 
much  simpler  picture  arises  and  consequently  more  accurate  fits  are  possible.  Conceptually,  one  could  also 
improve  the  fits  by  using  a  larger  number  of  parameters  in  the  model,  however,  this  would  defeat  the  purpose 
of  computational  efficiency  that  is  the  basis  for  seeking  reduced  kinetics. 

C.  Recovery  of  the  energetics 

According  to  eq.  8,  to  recover  the  value  of  T  in  the  reduction  scheme,  the  heavy  species  model  should  focus 
on  an  appropriate  representation  of  CP}h  and  I\h ■  Plots  (not  shown)  of  CPth  at  various  po  values,  calculated 
using  the  LLNL  model  over  a  wide  range  of  <j>,  show  that  these  curves  exhibit  a  very  modest  variation  over 
the  range  of  strong  Nc  decay;  moderate  (i.e.  up  to  50%)  variation  occurs  only  after  Nc  values  are  very 
small  to  negligible.  This  indicates  that  the  T  recovery  is  primarily  governed  by  A'/t  as  far  as  heavy  species 
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modeling  is  concerned. 

Examples  of  K/,  model  fits  to  the  LLNL  calculated  values  are  shown  in  Figs.  5  and  6.  First,  it  is  clear 
that  the  plots  resemble  functionally  those  of  Kc,  which  is  not  surprising  given  the  definitions  of  eqs.  3  and 
11.  Second,  the  energetics  is  remarkably  well  recovered  for  a  variety  of  (p  and  Ts  (Fig.  5),  including  very 
rich  mixtures  (Fig.  6).  Similarly,  excellent  to  very  good  results  are  obtained  for  lean  mixtures  (not  shown). 
Should  more  accurate  Kh  fits  be  desirable,  the  Kc  discussion  on  this  topic  is  entirely  pertinent. 

IV.  Completing  the  kinetic  model 

While  a  detailed  discussion  of  the  remaining  kinetic-model  parts  is  not  presented,  the  process  to  complete 
the  model  is  clear.  The  model  is  completed  with  the  fitting  of  quasi-steady  gain  rates  KGi  (where  i  stands  for 
H20,  C02,  02,  H,  CO,  H2,  CH4,  H202,  C2H2,  C2H4,  CH20),  of  quasi-steady  loss  rate  KLt  (where  i  stands 
for  H20,  02,  H,  H2,  this  rate  being  null  for  the  other  light  species)  and  of  the  quasi-steady  light  species  mole 
fraction  Xi.  To  show  that  the  Kc  fitting  has  set  a  methodology  for  completing  the  model,  displayed  in  Fig.  7 
is  the  quasi-steady  rate  function  KG  for  CO  at  po  =  30  bar.  The  plot  was  obtained  using  the  LLNL  kinetics. 
Function  KGco  represents  the  net  gain  of  molar  density  Nco  due  to  all  rates  involving  the  heavy  species 
group  (note:  KLco  is  null).  The  plot  graphically  resembles  Kc  and  Kh,  indicating  a  consistent  quasi-steady 
functional  behavior  for  the  heavy  group  contribution  similar  to  that  previously  identified  for  Kc.  Moreover, 
with  increasing  <p  beyond  (p  =  1,  an  asymptotic  self-similar  behavior  is  evident  in  the  fast  reaction  regime. 
This  plot  is  typical  of  all  KGi  rates  and  shows  that  the  same  fitting  methodology  used  for  Kc  and  Kf,  will 
apply. 

Regarding  the  rates  KLt  (except  for  possibly  i  =H),  they  also  have  the  specific  graphical  form  (not  shown) 
seen  for  Kc,Kh  and  KGco ;  thus,  for  constructing  KLj  and  X;h  fits,  the  same  basic  fitting  methodology  will 
apply  as  for  Kc  and  Kh- 


V.  Conclusions 

The  LLNL  heptane  kinetic  model  has  been  examined  focussing  on  a  potential  very  computationally 
efficient  kinetic  reduction  scheme.  Species  were  categorized  into  light  or  heavy  (C  number  n  ^  3).  The  heavy 
species  were  decomposed  into  defined  13  constituents,  and  it  was  shown  that  a  non-dimensional  temperature 
9  depending  on  the  pressure  variable  N*  and  the  equivalence  ratio  cp  can  be  defined  so  that  the  ratio  of  the 
total  constituent  molar  density  Nc  to  the  product  of  equivalence  ratio  and  a  nondimensional  pressure  N* 
exhibits  a  self-similar  behavior  over  a  wide  range  of  initial  po ,  <p  and  a  model-defined  temperature  Ts  entering 
the  computation  of  9.  The  number  of  progress  variables  is  thus  reduced  from  160  progress  variables  to  12 
(11  for  perfectly  stirred  reactor  cases  where  H  is  also  quasi-steady):  the  unsteady  light  species  and  Nc. 

Assessment  of  the  total  constituent  reaction  rate,  Kc .  behavior  revealed  that  it  is  quasi-steady  and 
features  two  distinct  regimes:  incubation  and  fast  reaction.  The  rate  was  fitted  in  terms  of  a  total  of  four 
parameters:  the  physical  parameters  T  or  9  (a  surrogate  for  T),  (p ,  N*  and  a  modeling  parameter  Ts ;  9  itself 
is  a  function  of  cp,  N*  and  Ts .  The  same  functions  were  used  to  achieve  the  fits  over  the  entire  region  covered 
by  the  four  parameters.  These  plots  reproduced  with  excellent  to  very  good  accuracy  the  reaction  rate  over 
the  significant  lifetime  of  the  total  constituent  and  for  a  range  of  po  and  cp  that  is  much  wider  than  the  typical 
one  used  for  displaying  reduced  mechanisms.  The  energetics  associated  with  the  heavy  species  reaction  rate, 
Kh,  was  also  recovered  through  fits  which  were  functionally  similar  to  those  for  Kc.  Expectably,  it  is  found 
that  the  Kj,  value  rather  than  the  value  of  the  heavy  species  heat  capacities  at  constant  pressure,  Cp 
contributes  mostly  to  the  enthalpy  change  associated  with  the  heavy  species  reaction.  The  CPth  also  requires 
fitting.  Mathematically,  the  shown  fits  achieved  substantial  accuracy  over  seven  rate  decades  considering  the 
fact  that  they  depend  upon  four  variables. 

The  model  is  completed  once  the  gain/loss  of  light  species  molar  density  from  the  heavy  species  and  the 
mole  fraction  Xi  (relative  to  the  light  group)  of  quasi-steady  light  species  are  fitted.  Examination  of  these 
gain/loss  rates  shows  that  they  also  have  a  quasi-steady  behavior.  As  an  example,  it  was  shown  that  the 
gain-rate  functional  form  with  respect  to  9  emulates  those  of  Kc  and  Kh,  making  the  modeling  strategy 
clear.  This  finding  highlights  the  fact  that  the  fitting  technique  provides  a  methodology  which  can  be 
repeatedly  used  to  obtain  an  accurate  representation  of  full  or  skeletal  kinetic  models.  The  advantage  of  the 
proposed  model  is  that  the  kinetic  reduction  is  performed  before  embedding  the  kinetics  in  a  turbulent  flow 
computation  and  would  be  used  as  a  library  of  rates  rather  than  needing  ‘on-the-fly’  kinetic  rate  calculations. 
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Summarizing,  the  proposed  model  features  12  ultimate  progress  variables  describing  heptane  oxidation. 
The  model  requires  fitting  16  quasi-steady  rate  functions,  11  curves  for  light  quasi-steady  mole  fractions,  and 
Cp,h-  Also  required  are  conventional  rates  between  members  of  the  light  group,  which  are  available  in  the 
literature.  Because  this  model  is  based  on  the  Nc  rate  rather  than  on  that  of  individual  heavy  species,  even  if 
the  number  of  species  increases  with  increased  C  number  in  the  alkane  group,  providing  that  the  quasi-steady 
rate  aspect  persists,  then  extension  of  this  model  to  higher  alkanes  should  be  conceptually  straightforward; 
although  it  remains  to  be  seen  if  the  functional  fits  would  remain  valid  or  would  require  reconstruction. 
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Figure  2.  Kc  versus  6  for  po  =  20  bar.  Symbols  □,  >,A,0>V,<  represent  selected  data  from  the  LLNL  runs; 

lines  •••,  —  •  —  •  — -  -, - , - —  represent  the  corresponding  fits.  Temperatures  units  are  in 

degrees  K.  The  legend  is:  •  ••</>  =  0.5,  Ts  =  825; - <p  —  0.5,  Ts  =  715,  —  •  —  •  —  (p  =  1.0,  Ts  =  825; - <p  =  1.0, 

Ts  =  715; - <p  =  2.0,  Ts  =  825;  -  ••-••-0  =  2.0,  Ts  =  715. 


Figure  3.  Kc  versus  6  for  (p  —  1/3  at  different  po  and  Ts.  Symbols  □,  >,A,0?V,<]  represent  selected  data  from 

the  LLNL  runs;  lines  - , - , - ,  —  —  represent  the  corresponding  fits.  Units:  pressures 

in  bar,  temperatures  in  degrees  K.  The  legend  is:  •  •  •  po  =  10,  Ts  =  825; - po  =  10,  Ts  =  715,  —  •  —  •  —  po  =  20, 

Ts  =  825; - po  =  20,  Ts  =  715; - po  =  40,  Ts  =  825;  p0  =  40,  Ts  =  715. 


Figure  4.  Kc  versus  6  for  a  variety  of  conditions.  Symbols  □,  >,A,05^75<1  represent  selected  data  from  the 

LLNL  runs;  lines  •••,  —  •  —  •  — , - , - , - ,  —  —  represent  the  corresponding  fits.  Units:  pressures 

in  bar,  temperatures  in  degrees  K.  (a)  <p  =  4,  (•  •  •  po  =  10,  Ts  =  825; - po  =  10,  Ts  =  715,  —  •  —  •  —  po  =  20, 

Ts  =  825; - p0  =  20,  Ts  =  715; - p0  =  40,  Ts  =  825;  Po  =  40,  Ts  =  715).  (b)  p0  =  10  and  <p  =  1  (•  •  • 

Ts  =  825; - Ts  =  715,  -  •  -  •  -  Ts  =  655; - Ts  =  560). 
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Figure  5.  K ^  versus  6  at  fixed  po  =  20  bar  for  different  (ft  and  Ts  values  corresponding  to  (a)  Fig.  2a,  and  (b) 
Fig.  2b.  The  legend  is  in  Fig.  2  caption. 


Figure  6.  K ^  versus  9  for  (a)  (ft  =  4  corresponding  to  Fig.  4a,  and  (b)  (ft  =  1  and  po  =  10  bar  corresponding  to 
Fig.  4b.  The  legend  is  in  Fig.  4  caption. 


0 

Figure  7.  KG  for  CO  versus  6.  po  is  in  bar.  Plots  obtained  using  the  LLNL  data  in  CHEMKIN  II. 
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Computationally  efficient  modeling  of  n-heptane  oxidation 
using  constituents  and  species 
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Abstract 

A  simplified  model  is  proposed  for  the  kinetics  of  alkane  oxidation  in  air,  based  on  a  decomposition  of  heavy 
(carbon  number  greater  or  equal  to  3)  hydrocarbons  into  a  13  constituent  radical  base.  The  behavior  of  this  base 
is  examined  in  test  computations  for  n-heptane  utilizing  Chemkin  II  with  LLNL  data  inputs,  placing  emphasis  on 
modeling  to  predict  the  heat  release  and  temperature  evolution.  A  normalized  temperature  was  constructed  which 
when  used  to  plot  the  total  constituent  molar  density  divided  by  the  product  of  the  equivalence  ratio  and  a  nondi- 
mensional  pressure,  reveals  a  self-similar  behavior  of  the  plotted  variable  over  a  wide  range  of  initial  pressures  and 
equivalence  ratios.  Examination  of  the  LLNL  kinetics  shows  that  the  total  constituent  molar  density  rate  follows 
a  quasi-steady  behavior.  This  reaction  rate  was  curve  fitted  along  with  the  corresponding  enthalpy  production. 
The  fits  are  shown  against  the  normalized  temperature  for  various  equivalence  ratios  and  initial  nondimensional 
pressures  and  comparisons  with  the  LLNL  kinetics  are  very  favorable.  The  model  reduces  the  LLNL  n-heptane 
mechanism  from  160  species  (progress  variables)  and  1540  reactions  to  12  progress  variables,  16  quasi-steady 
rates  (associated  with  heavy  species),  162  conventional  reaction  rates  (light  species)  and  1 1  other  functional  forms 
(i.e.  fits  for  the  mean  heavy-species  heat  capacity  at  constant  pressure,  the  enthalpy  release  rate  of  the  heavy 
species,  and  the  molar  fraction  of  quasi-steady  light  species).  The  proposed  kinetic  mechanism  is  valid  over  a 
pressure  range  from  atmospheric  to  60  bar,  temperatures  from  600  K  to  2500  K  and  equivalence  ratios  from  0.125 
to  8.  This  range  encompasses  diesel,  HCCI  and  gas  turbine  engines,  including  cold  ignition;  and  NO;,,,  CO  and 
soot  pollutant  formation  in  the  lean  and  rich  regimes,  respectively. 

Keywords:  reduced  n-heptane  kinetics 


1.  Introduction 

Computational  models  for  turbulent  reactive  flows 
present  challenges  that  have  yet  to  be  frilly  met. 
First,  combustion  of  typical  hydrocarbons  proceeds 
through  the  formation  of  hundreds  of  species  interact¬ 


ing  through  thousands  of  reactions.  Second,  a  com¬ 
plete  numerical  resolution  of  high  Reynolds  number 
flow  requires  a  very  large  number  of  grid  points  that 
is  beyond  current  typical  supercomputer  capabilities. 
Thus,  combining  detailed  reaction  rate  models  with 
the  fundamental  flow  equations  yields  a  numerical 


problem  that  is  computationally  prohibitive.  A  dras¬ 
tic  (at  least  order  of  magnitude)  reduction  in  the  di¬ 
mension  of  reaction  kinetics  is  required  for  a  feasible 
computational  model.  In  spirit,  this  reduction  is  simi¬ 
lar  to  that  required  in  Large  Eddy  Simulations  (where 
energetically  significant  flow  scales  are  computed  and 
the  others  are  modeled)  with  respect  to  Direct  Numer¬ 
ical  Simulations  (where  all  flow  scales  are  computed). 
Since  each  species  has  a  characteristic  time  scale,  in 
the  kinetic  reduction  it  is  also  desirable  to  compute 
only  those  species  having  an  essential  characteristic 
time  scale  (to  be  defined)  and  model  the  kinetics  of 
the  remaining  species. 

We  categorize  all  species  involved  in  hydrocar¬ 
bon  oxidation  into  light  species  and  heavy  species; 
the  heavy  species  are  those  having  a  carbon  number 
n  ^  3.  Based  on  this  definition,  a  method  for  simpli¬ 
fying  alkane  combustion  is  proposed  and  tested  for 
n-heptane.  The  method  focusses  on  the  total  con¬ 
stituent  radical  molar  density  (see  Section  2  defin¬ 
ition  and  constituent  radical  base  formation)  of  the 
heavy  species  rather  than  on  the  much  larger  num¬ 
ber  of  these  species.  That  is,  rather  than  following  all 
species  through  their  reaction  coordinates,  we  follow 
a  reduced  set  of  reaction  coordinates;  this  reduced  set 
is  called  a  base.  Our  primary  objective  is  the  accurate 
determination  of  the  energetics,  i.e.  heat  release  and 
temperature  evolution;  at  this  initial  stage  of  model 
development,  neither  NCL  nor  soot  kinetics  are  ad¬ 
dressed.  In  Section  3,  the  behavior  of  the  base  set  is 
examined  in  the  context  of  a  detailed  reaction  model 
for  heptane  combustion  as  given  by  LLNL  [1]  and  the 
proposed  model  is  tested  against  the  LLNL  kinetic 
scheme.  Section  4  is  devoted  to  a  discussion  on  the 
completion  of  the  model.  Finally,  in  Section  5  a  sum¬ 
mary  of  the  accomplishments  and  a  perspective  on  the 
future  work  are  offered. 

2.  Species  representation:  a  reduced  base  set 

2.1.  Heavy  species  decomposition  into  constituent 
radicals 

Consider  the  count  of  atoms  (H,C,0,N)  in  species 
i.  We  define  Wji  as  the  number  of  species  j  in  species 
i  that  produces  the  correct  atom  count,  where  the  in¬ 
dex  j  denotes  the  species  set  of  air  (02,N2)  and  final 
combustion  products  (H20,CC>2).  (This  means  that 
mathematically,  we  consider  species  i  as  being  com¬ 
posed  of  air  and  final  combustion  product  species.) 
The  molar  densities  of  species  in  this  set  are 

Nj=N°  -J^  wjiNi  (1) 

where  Nj  is  the  effective  total  molar  density  of 
species  j.  The  same  exercise  of  atom  counting  and 
species  representation  by  the  j  set  can  be  performed 
for  all  species  in  the  mixture  to  obtain  Nj  .  That  is, 
the  atomic  elements  molar  densities  are  represented 
by  species  j  molar  densities.  Conservation  of  atoms 


in  the  reactions  makes  unnecessary  utilizing  rate  ex¬ 
pressions  for  reaction  coordinates  of  air  species  and 
complete-combustion  products,  as  a  one-to-one  map 
exists  between  these  and  the  atom  set.  The  heavy 
species  constituent  radicals  (CH2,  CHs,  CH,  C2H3, 
C2H2,  C2,  HC2,  CO  (keto),  HCO,  HO,  H02,  00,  O) 
form  a  set  of  1 3  entities  from  which  any  heavy  species 
or  radical  may  be  constructed.  The  total  constituent 
molar  density  is 
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NC  =  J2  Nk  (2) 

k=i 

where  the  molar  density  of  constituent  k  is  Nk  .  The 
decomposition  into  constituents  is  similar  to  group 
additivity  [2,  3],  except  that  it  considers  only  first  or¬ 
der  (compositional)  effects  and  only  interactions  with 
adjacent  groups.  The  reaction  rate  associated  with  Nc 
is  defined  by 


K  (3) 
at. 

We  define  a  reference  density  for  nitrogen  N2  from 
the  dry  air  value  at  pressure  pref  =  1  bar  and  Tre/= 
298.15  K,  giving  Nref—  31.5  mol/m3.  Using  Nref,  a 
dimensionless  N2  molar  density  is  defined  as 

N*  =  NN2/Nref  (4) 

which  acts  as  a  convenient  surrogate  variable  for  the 
pressure  p  as  the  N*  value  is  essentially  rate  invariant. 
That  is,  the  partial  pressure  of  N2  is  the  overwhelming 
contribution  to  p.  Basically,  the  ratio  of  N*  to  po  is 
nearly  constant.  A  normalized  temperature  variable  9 
is  defined  by 


6 

T-Ts 

Tr^,N *) 

(5) 

Tr 

ee  2065(7V*)°'0M</>) 

(6) 

w{4>) 

,  1.5  +  1.310 

(V) 

1  +  0.710  -f 

where  T„  is  a  modeling  parameter  representing  the 
lowest  temperature  at  which  I\c  >  0.  The  6  defi¬ 
nition  is  such  that  Nc  decreases  during  stoichiomet¬ 
ric  combustion  from  its  initial  value  by  three  orders 
of  magnitude  (delving  into  higher  order  of  magnitude 
decrease  runs  the  risk  of  encountering  round-off  and 
truncation  errors)  for  9  >  0.6,  a  value  chosen  so  that 
all  9  values  remain  below  unity  for  all  test  case  calcu¬ 
lations. 

The  contribution  of  the  heavy  group  to  the  rate  of 
change  of  a  light  species  molar  density  Nz  with  mole 
fraction  Xt  (relative  to  the  light  group)  is  expressed 
in  terns  of  quasi-steady  gain  and  loss  rates  KGi  and 
K L,  as 


dNi 


Nc(KGi  -  XiKLi). 


(8) 


Aside  from  the  molar  reaction  rates,  determining  the 
temperature  evolution  in  a  reactive  system  requires 
knowledge  of  the  species  molar  enthalpies  and  heat 
capacities.  For  species  i,  the  molar  enthalpy  may 
be  expressed  as  hi  =  h®  +  hfT)  where  h°  is 
the  heat  of  formation,  hi  is  the  sensible  enthalpy 
and  T  is  the  temperature.  The  heat  of  formation 
is  taken  at  the  above  reference  conditions,  giving 
hi  =  fj  Cv  idT  where  Cv  ;  is  the  molar  constant- 
pressure  heat  capacity.  A  heat  of  combustion  for 
species  i  is  given  by 


2.3.  Final  base  set 

The  reaction  base  set  is  composed  of  a  global  (i.e. 
total)  constituent  radical  and  1 1  light  molecules  or 
free  radicals  of  the  first  subset  (Table  1).  Thus,  there 
are  only  12  progress  variables.  To  compute  Nc  one 
must  model  Kc ;  to  compute  the  1 1  light  species,  one 
must  model  the  KGi ,  the  ALL;,  the  X,  of  the  sec¬ 
ond  subset  and  the  conventional  light  species  reaction 
rates  of  the  first  subset. 

3.  Results 


hi  =  h°  u\;ihj  (9) 

3 

under  the  assumption  that  water  remains  in  vapor 
state.  The  energy  equation  is 


NcCp,h+  ]T  Cp.iNi 

i£lights 

—  hiTZi  +  Nc(RuTref)Kh  (10) 

i£lights 

where 

7 Zi  =  (dNi/dt)reac,  (11) 

Ru  is  the  universal  gas  constant, 

Cp.h  =  (  ]T  Cp,iNi)/Nc,  (12) 

l£heavies 

and  the  rate  of  enthalpy  change  due  to  heavy  species 
l  chemical  kinetic  changes  is  defined  as 


Kh  =  - 


1 

RuTrefNc 


(13) 


2.2.  Light  species  as  part  of  the  base 

Light  species  are  not  subject  to  meaningful  decom¬ 
position  according  to  the  above  procedure.  Examina¬ 
tion  of  runs  made  using  the  LLNL  database  indicate 
that  the  light  species  should  be  categorized  in  two 
subsets.  The  species  in  the  first  subset  require  rate 
equations  and  are  part  of  the  base  set.  The  species  in 
the  second  subset  have  a  quasi-steady  behavior  and 
require  fits  of  their  mole  fraction  value  as  a  func¬ 
tion  of  the  state  variables  (T,  N* ,  <j>)  and  modeling 
parameters  (here,  only  Ts).  This  first  subset  consists 
of:  H20,  C02,  02,  H,  CO,  H2,  CH4,  H202,  C2H2, 
C2H4,  CH20.  The  second  subset  are  the  radicals:  O, 
CH,  CH2,  CH3,  HO,  HCO,  H02,  HC2,  C2H3.  Thus 
this  entire  set  consists  of  20  light  species.  For  their 
energetics,  values  of  hc  are  fitted  from  the  literature 
and  Cp’ s  are  obtained  as 

Cp/Ru  =ah  +  bh  In (T/Tref).  (14) 

Other  properties  of  the  light  molecules  and  light  free¬ 
form  of  the  constituent  radicals  are  given  in  Table  1 
[4-7]. 


3.1.  Identification  of  a  self-similar  behavior 

Showing  the  appropriateness  of  the  scaling  with  9 

rather  than  using  T,  the  evolution  of  Nc  scaled  by 
(0  x  N*)  for  fixed  po  at  several  0  values  and  at  fixed 
<f>  for  several  po  values  is  illustrated  in  Fig.  1  using 
the  LLNL  kinetics.  The  normalization  attainable  with 
6  scaling  is  clear  and  so  is  the  fact  that  7Vc/(0  x  N*) 
is  a  similarity  variable.  At  fixed  p o  =  20  bar,  Fig.  la, 
all  Nc/(<j>  x  N*)  curves  coincide  nearly  on  the  same 
curve,  with  the  somewhat  dissenting  curve  for  0  =  4 
showing  that  for  very  rich  mixtures  there  are  indeed 
constituents  left  after  the  oxidizer  was  depleted.  How¬ 
ever,  even  for  as  rich  a  mixture  as  0  =  2,  similar¬ 
ity  holds,  making  this  model  valid  in  realms  beyond 
those  in  advanced  reduced  schemes  [8]  where  the  up¬ 
per  limit  of  the  scheme  validity  is  0  =  1 . 5,  at  most  [9] 
0  =  2.0  or  exceptionally  [10]  0  =  3.  The  behavior 
in  Fig.  la  is  typical  of  all  fixed  po  plots.  In  Fig.  lb, 
corresponding  plots  at  stoichiometric  conditions  are 
shown,  which  are  generic  of  all  fixed  0  plots.  Again, 
the  similarity  achieved  with  the  two  normalized  vari¬ 
ables  is  remarkable,  and  only  at  a  combination  of  high 
po  values  and  small  Ta  values  (typical  of  the  lowest  T 
limit  at  which  Kc  rates  are  finite),  are  small  to  mod¬ 
erate  departures  from  the  excellent  normalization  per¬ 
ceptible. 

We  note  that  by  6  ~  [0.55, 0.6],  Nc  ~  0,  meaning 
that  modeling  errors  past  that  6  range  are  somewhat 
irrelevant  to  the  model’s  accuracy  since  the  reaction 
is  practically  finished.  Although  it  would  be  tempt¬ 
ing  to  curve  fit  dNc/dt  from  Fig.  1  and  adopt  it  as  a 
representative  kinetic  rate,  the  fact  is  that  the  scaling 
hides  the  precise  detailed  formulation  that  is  required 
for  a  sufficiently  accurate  kinetics  to  predict  T.  Thus, 
fits  for  Kc  are  instead  needed  to  achieve  appropriate 
accuracy,  as  discussed  below. 

3.2.  Global  constituent  molar  density  model 

Plots  of  Kc  versus  6  show  that  there  are  two  differ¬ 
ent  regions  of  the  curve.  For  small  9  values,  there  is 
a  combustion  incubation  region,  with  reaction  times 
values  in  ms  units.  As  9  increases,  one  encounters 
the  fast  reaction  region  with  reaction  times  values  in 
ps  units.  Numerous  tests  made  with  the  LLNL  ki¬ 
netics  showed  that  Kc  represents  only  a  small  dif¬ 
ference  between  production  and  loss  rates  which  are 
nearly  equal  (not  shown)  and  are  much  larger  than 
Kc.  The  Kc  value  is  typically  about  1.5%  to  5%  of 


either  production  or  loss  rate.  The  maximum  differ¬ 
ence  between  production  and  loss  rates  is  about  5% 
for  N*  <  1,  then  decreases  to  about  3%  at  larger 
N* .  These  maximum  differences  tend  to  occur  dur¬ 
ing  very  large  values  of  dT / dt.  This  fact  suggests  that 
Kc  could  be  treated  as  a  rate  of  a  nearly  quasi-steady 
process.  If  fits  of  Kc  can  be  obtained  over  a  signifi¬ 
cant  region  of  parametric  behavior,  this  model  would 
be  compact  and  represent  a  considerable  reduction  in 
the  problem  dimensionality.  Noteworthy,  had  the  ap¬ 
proach  been  to  find  Kc  from  slopes  of  Nc  in  Fig.  1, 
it  is  clear  that  the  important  incubation  region  would 
have  been  missed. 

Fitting  Kc  as  a  function  of  T  and  parameters  N* , 
0  and  Ta  presented  a  major  challenge  because  of 
the  non-monotonic  Kc  behavior  and  also  because  Kc 
spans  so  many  decades.  Fits  were  obtained  by  first 
inspecting  the  LLNL  kinetics  to  understand  the  rate 
constant,  K ,  variation.  This  variation  showed  an  in¬ 
cubation  region  ( 8  <  0.2)  featuring  first  a  maximum, 
Kmx ,  followed  by  a  minimum,  Kmn ,  and  a  high-T 
region  (6  >  0.2)  where  Kc  monotonically  increases. 
Thus,  the  functional  fits  are  performed  in  three  sep¬ 
arate  regions,  with  connection  constraints  between 
consecutive  regions  (the  ultimate  fits  are  piecewise 
continuous).  The  non-monotonic  incubation- region 
behavior  is  captured  by  using  a  cubic  functional  map¬ 
ping  from  K(T)  to  y(T )  through  the  form 

K(T)  =  i( Kmx  +  Kmn) 

+  -^(Kmx  —  Kmn)(y3  —  3y)  (15) 

where  y(T)  =  —  1  at  T  =  Tmx  which  is  the  loca¬ 
tion  of  Kmx  and  y(T)  =  1  at  T  =  Tmn  which  is 
the  location  of  Kmn  ■  In  fact,  the  range  of  y  is  from 
slightly  smaller  than  -2  at  T  =  Ta  (i.e.  Kc  >  0)  to 
y  <  1.8.  The  values  of  Kmx,Kmn,Tmx  and 

Tmn 

are  fitted  as  polynomials  in  parameters  ln(7V*),  ln(0) 
and  ln(Ts/Tre/).  Also  fitted  is  To  which  is  the  T 
value  at  which  y  =  0.  Two  different  changes  of 
variables  are  made  from  y(T)  to  y(z)  depending  on 
how  T  compares  to  To.  For  T  <  To,  a  variable 
z  =  (To  —  T)/(To  —  Tmx)  is  defined,  which  means 
that  y  =  —  1  at  z  =  1.  To  achieve  the  mapping,  we 
first  generate  (y,  z)  pairs  which  are  functions  of  pa¬ 
rameters  ln(lV*),  ln(0)  and  ln(Ts/Tre/).  To  match 
these  pairs,  a  ratio  of  appropriate  functions  (e.g.  poly¬ 
nomials  a  +  bz  +  cz2  +  ...;  power  functions  of  type 
az;  and  combinations  of  the  two  functional  forms)  is 
chosen  to  produce  the  y(z)  mapping.  For  T  >  To, 
i.e.  y  ^  0,  a  similar  procedure  is  used  for  y(z)  where 
now  z  =  (T  —  To)/(Tm„  —  To).  Beyond  incuba¬ 
tion,  for  the  high-T  region,  ln(Lf)  is  fitted  as  a  fifth 
order  polynomial  in  6P  where  p  =  1  for  0  ^  0.5  and 
p  =  (0.86  +  O.280)/(O.74  +  0.520)  for  0  >  0.5. 
Scrutiny  of  the  LLNL  kinetics  showed  that  the  largest 
value  of  T,  is  825  K.  Then,  for  Ts  =  825  K,  the 
polynomial  coefficients  were  fitted,  separately  for  the 
lean/stoichiometric  (0  ^  1)  and  rich  (0  >  1)  regions, 


as  polynomials  in  ln(iV*)  and  ln(0).  To  obtain  fits 
for  Ts  <  825  K,  a  scaling  procedure  was  used  from 
the  Ts  =  825  K  fits  to  match  particular  K  values  at  a 
specific  point  in  the  incubation  region. 

All  fits  further  presented  in  Figs.  2-6  have  been 
achieved  with  the  fixed  set  of  mapping  functions  de¬ 
scribed  above.  That  is,  no  local  corrections  of  any  sort 
were  used,  nor  exceptions  from  those  functions  were 
considered.  Local  corrections  in  a  four-dimensional 
space  are  not  computationally  practical. 

Curve  fits  to  Kc  computed  using  the  LLNL  kinet¬ 
ics  are  plotted  in  Figs.  2-4.  In  Fig.  2a  and  2b  (fixed 
po,  various  0  and  Ts),  and  Fig.  3a  and  3b  (fixed  0  in 
the  lean  regime,  various  p o  and  Ts)  we  differentiate 
between  the  incubation  and  fast-reaction  regions  dis¬ 
cussed  above  and  show  on  each  figure  up  to  5  decades 
of  kinetic  rates.  On  each  figure,  symbols  represent 
the  LLNL-based  computation  and  the  symbol  loca¬ 
tion  represents  a  selected  sample  of  time  stations  dur¬ 
ing  the  calculation.  The  curve  fits  are  shown  as  lines 
and  are  based  on  time  stations  at  fixed  intervals  of  Ad. 
The  agreement  on  the  logarithmic  scale  is  excellent  to 
very  good  and  deterioration  of  the  fits  at  9  >  0.55  cor¬ 
respond  to  very  small  values  of  Nc,  as  noted  above. 
Particularly  important  is  the  excellent  agreement  be¬ 
tween  fits  and  LLNL  kinetics  in  the  region  off?  <  0.2 
representing  54%  of  the  Nc  decay  (see  Fig.l).  For 
9  >  0.2  the  agreement  is  still  excellent  as  long  as  Ts 
is  not  too  low;  good  agreement  still  occurs  in  the  sig¬ 
nificant  Nc  regions  (i.e.  8  <  [0.55,  0.6])  and  over  a 
much  wider  0  range  than  usually  shown  in  the  litera¬ 
ture.  Figure  4a  focusses  on  the  rich  regime,  showing 
the  entire  range  of  relevant  8  (see  Fig.  la)  at  vari¬ 
ous  po  and  T, .  The  agreement  of  fits  with  the  LLNL 
data  is  excellent,  bodying  well  for  future  studies  ad¬ 
dressing  soot  and  CO  formation.  Noteworthy,  the 
fits  encompass  a  wide  p o  range  including  the  desired 
high-pressure  situations  for  engine  combustion,  and  a 
much  wider  0  range  than  typically  used  in  kinetic  re¬ 
duction  schemes  [8]  where  generally  0  6  [0.5, 1.5], 
sometimes  [9]  0  €  [0.5,  2.0]  and  exceptionally  [10] 
0  £  [0.7, 3.0].  To  show  the  influence  of  T„,  in  Fig. 
4b  is  depicted  the  case  for  which  the  widest  range  of 
Ts  was  available,  corresponding  to  po  =  10  bar  and 
0=1.  Similar  to  the  other  figures,  the  agreement 
with  the  LLNL  kinetics  is  excellent  to  very  good  up 
to  9  <  0.2  and  somewhat  deteriorates  past  that  sta¬ 
tion. 

This  deterioration  is  a  manifestation  of  the  need  for 
matching  two  different  functional  forms:  one  in  the 
incubation  region  (y  >  0),  and  the  other  in  the  high- 
T  region.  The  matching  point  is  at  j/max;  as  2/max 
increases,  the  high-T  region  fit  improves  while  the 
incubation  fit  worsens,  and  vice  versa  for  decreasing 
?/max.  Fits  can  also  be  improved  if  one  is  less  ambi¬ 
tious  in  covering  the  large  0  range  treated  here.  The 
functional  behavior  of  0  outliers  is  somewhat  differ¬ 
ent  from  that  for  moderate  0.  For  0  <  1  and  small,  the 
high-T  rates  are  quite  modest  and  below  Kmx.  For 
0  >  1,  as  0  increases,  the  ratio  KmX/Kmn  decreases 
and  high-T  rates  exceed  Kmx  except  for  the  largest  0 


where  the  reaction  terminates  before  the  KmX  value 
is  reached  again  in  the  high-T  region.  This  behavior 
is  independent  of  po  .  This  is  to  say  that  if  the  <f>  out¬ 
liers  are  eliminated,  a  much  simpler  picture  arises  and 
consequently  more  accurate  fits  are  possible.  Con¬ 
ceptually,  one  could  also  improve  the  fits  by  using  a 
larger  number  of  parameters  in  the  model,  however, 
this  would  defeat  the  purpose  of  computational  effi¬ 
ciency  that  is  the  basis  for  seeking  reduced  kinetics. 

3.3.  Recovery  of  the  energetics 

According  to  eq.  10,  to  recover  the  value  of  T  in 
the  reduction  scheme,  the  heavy  species  model  should 
focus  on  an  appropriate  representation  of  Cv,h  and 
Kh.  Plots  (not  shown)  of  Cp,h  at  various  po  values, 
calculated  using  the  LLNL  model  over  a  wide  range 
of  <j>,  show  that  these  curves  exhibit  a  very  modest 
variation  over  the  range  of  strong  Nc  decay;  moderate 
(i.e.  up  to  50%)  variation  occurs  only  after  Nc  values 
are  very  small  to  negligible.  This  indicates  that  the  T 
recovery  is  primarily  governed  by  Kh  as  far  as  heavy 
species  modeling  is  concerned. 

Examples  of  Kh  model  fits  to  the  LLNL  calculated 
values  are  shown  in  Figs.  5  and  6.  First,  it  is  clear  that 
the  plots  resemble  functionally  those  of  Kc,  which  is 
not  surprising  given  the  definitions  of  eqs.  3  and  13. 
Second,  the  energetics  is  remarkably  well  recovered 
for  a  variety  of  <j>  and  Ta  (Fig.  5),  including  very  rich 
mixtures  (Fig.  6).  Similarly,  excellent  to  very  good 
results  are  obtained  for  lean  mixtures  (not  shown). 
Should  more  accurate  Kh  fits  be  desirable,  the  Kc 
discussion  on  this  topic  is  entirely  pertinent. 

4.  Completing  the  kinetic  model 

Brevity  does  not  permit  the  detailed  discussion  of 
the  remaining  kinetic -model  parts,  but  the  process  to 
complete  the  model  is  clear.  What  remains  to  be  done 
is  the  fitting  of  quasi-steady  gain  rates  KGi  (where 
i  stands  for  H20,  C02,  02,  H,  CO,  H2>  CH4,  H202, 
C2H2,  C2Fl4,  CH20),  of  quasi-steady  loss  rate  KLi 
(where  i  stands  for  H20,  02,  H,  H2 ,  this  rate  be¬ 
ing  null  for  the  other  light  species)  and  of  the  quasi¬ 
steady  light  species  mole  fraction  X;.  To  show  that 
the  Kc  fitting  has  set  a  methodology  for  completing 
the  model,  displayed  in  Fig.  7  is  the  quasi-steady  rate 
function  KG  for  CO  at  po  =  30  bar.  The  plot  was 
obtained  using  the  LLNL  kinetics.  Function  KGco 
represents  the  net  gain  of  molar  density  Nco  due 
to  all  rates  involving  the  heavy  species  group  (note: 
KLco  is  null).  The  plot  graphically  resembles  Kc 
and  Kf,  ,  indicating  a  consistent  quasi-steady  func¬ 
tional  behavior  for  the  heavy  group  contribution  sim¬ 
ilar  to  that  previously  identified  for  K,,.  Moreover, 
with  increasing  <f>  beyond  <f>  =  1,  an  asymptotic  self¬ 
similar  behavior  is  evident  in  the  fast  reaction  regime. 
This  plot  is  typical  of  all  KGi  rates  and  shows  that  the 
same  fitting  methodology  used  for  Kc  and  I\h  will 
apply. 

Regarding  the  rates  KLi  (except  for  possibly  i  = 
FI),  they  also  have  the  specific  graphical  form  (not 


shown)  seen  for  Kc,Kh  and  KGco ;  thus,  for  con¬ 
structing  KLi  and  X,  fits,  the  same  basic  fitting 
methodology  will  apply  as  for  Kc  and  Kh  ■ 

5.  Conclusions 

The  LLNL  heptane  kinetic  model  has  been  exam¬ 
ined  focussing  on  a  potential  very  computationally  ef¬ 
ficient  kinetic  reduction  scheme.  Species  were  cate¬ 
gorized  into  light  or  heavy  (C  number  n  ^  3).  The 
heavy  species  were  decomposed  into  defined  1 3  con¬ 
stituents,  and  it  was  shown  that  a  non-dimensional 
temperature  8  depending  on  the  pressure  variable  N* 
and  the  equivalence  ratio  (f>  can  be  defined  so  that  the 
ratio  of  the  total  constituent  molar  density  Nc  to  the 
product  of  equivalence  ratio  and  a  nondimensional 
pressure  N*  exhibits  a  self-similar  behavior  over  a 
wide  range  of  initial  po,  <j>  and  a  model-defined  tem¬ 
perature  Ta  entering  the  computation  of  8.  The  num¬ 
ber  of  progress  variables  is  thus  reduced  from  160 
progress  variables  to  12(11  for  perfectly  stirred  reac¬ 
tor  cases  where  H  is  also  quasi-steady):  the  unsteady 
light  species  and  Nc. 

Assessment  of  the  total  constituent  reaction  rate, 
Kc,  behavior  revealed  that  it  is  quasi-steady  and  fea¬ 
tures  two  distinct  regimes:  incubation  and  fast  reac¬ 
tion.  The  rate  was  fitted  in  terns  of  a  total  of  four  pa¬ 
rameters:  the  physical  parameters  T  or  6  (a  surrogate 
for  T),  tp,  N*  and  a  modeling  parameter  Ta;  6  itself  is 
a  function  of  <j>,  N*  and  Ta.  The  same  functions  were 
used  to  achieve  the  fits  over  the  entire  region  covered 
by  the  four  parameters.  These  plots  reproduced  with 
excellent  to  very  good  accuracy  the  reaction  rate  over 
the  significant  lifetime  of  the  total  constituent  and  for 
a  range  of  po  and  <p  that  is  much  wider  than  the  typi¬ 
cal  one  used  for  displaying  reduced  mechanisms.  The 
energetics  associated  with  the  heavy  species  reaction 
rate,  Kh ,  was  also  recovered  through  fits  which  were 
functionally  similar  to  those  for  Kc.  Expectably,  it 
is  found  that  the  Kh  value  rather  than  the  value  of 
the  heavy  species  heat  capacities  at  constant  pressure, 
Cp,h,  contributes  mostly  to  the  enthalpy  change  as¬ 
sociated  with  the  heavy  species  reaction.  The  Cp,h 
also  requires  fitting.  Mathematically,  the  shown  fits 
achieved  substantial  accuracy  over  seven  rate  decades 
considering  the  fact  that  they  depend  upon  four  vari¬ 
ables. 

The  model  is  completed  once  the  gain/loss  of  light 
species  molar  density  from  the  heavy  species  and 
the  mole  fraction  X;  (relative  to  the  light  group)  of 
quasi-steady  light  species  are  fitted.  Examination  of 
these  gain/loss  rates  shows  that  they  also  have  a  quasi¬ 
steady  behavior.  As  an  example,  it  was  shown  that  the 
gain-rate  functional  form  with  respect  to  8  emulates 
those  of  Kc  and  I\h,  making  the  modeling  strategy 
clear.  This  finding  highlights  the  fact  that  the  fitting 
technique  provides  a  methodology  which  can  be  re¬ 
peatedly  used  to  obtain  an  accurate  representation  of 
foil  or  skeletal  kinetic  models.  The  advantage  of  the 
proposed  model  is  that  the  kinetic  reduction  is  per¬ 
formed  before  embedding  the  kinetics  in  a  turbulent 
flow  computation  and  would  be  used  as  a  library  of 


rates  rather  than  needing  ‘on-the-fly’  kinetic  rate  cal¬ 
culations. 

To  summarize,  the  proposed  reduced  model  fea¬ 
tures  12  ultimate  progress  variables  describing  hep¬ 
tane  oxidation.  The  reduced  model  requires  fitting 
1 6  quasi-steady  rate  functions,  1 1  curves  for  light 
quasi-steady  mole  fractions,  and  CPth-  Also  required 
are  conventional  rates  between  members  of  the  light 
group,  which  are  available  in  the  literature.  Because 
this  model  is  based  on  the  Nc  rate  rather  than  on 
that  of  individual  heavy  species,  even  if  the  number 
of  species  increases  with  increased  C  number  in  the 
alkane  group,  providing  that  the  quasi-steady  rate  as¬ 
pect  persists,  then  extension  of  this  model  to  higher 
alkanes  should  be  conceptually  straightforward;  al¬ 
though  it  remains  to  be  seen  if  the  functional  fits 
would  remain  valid  or  would  require  reconstruction. 
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bh 

h2 

0.0 

241.5 
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0.0 
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0.5663 

n2 

0.0 

0.0 

3.388 

0.469 
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mi 
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0.0 
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-241.5 
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co2 

-393.5 

0.0 

4.690 

1.390 
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473 

473 

2.50 

0.0 
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218.0 

339 

2.50 

0.0 

HO 

38 

159 

3.385 

0.3637 

HOO 

10.5 

131 

4.150 

1.307 

O 

249.2 

249.2 

2.536 

0.0 

CO 

-110.5 

283 

3.426 

0.4749 

HCO 

43.1 

558 

4.154 

1.2875 

ch4 

-74.6 

802 

3.797 

4.305 

ch3 

146 

902 

4.440 

2.249 

ch2 

390 

1025 

3.973 

1.3015 

CH 

596 

1110 

3.220 

0.7136 

c2h3 

300 

1449 

5.1 

3.5 

hc2 

566 

1474 

4.434 

1.404 

c2 

838* 

1625 

4.58 

0.0 

NO 

90.3 

90.3 

3.533 

0.4508 

no2 

33.2 

33.2 

4.691 

1.151 

h2o2 

-136 

106 

5.269 

1.880 

HCOH 

-109 

526.5 

4.27 

2.546 

c2h2 

228 

1257 

5.368 

2.294 

c2h4 

52.5 

1323 

5.383 

4.676 

*  Alternate  value  of  832  from  the  CRC  tables. 


Table  1 :  Thermodynamic  properties  of  molecules  and  free 
radicals.  h°  and  hc  (heats  of  formation  and  combustion,  re¬ 
spectively),  in  kJ/mol;  constants  ah  and  bh  for  molar  heat 
capacity  in  the  form  Cv /Ru  =  ah  -\-bhln(T/Tref);Tref 
=  298.15  K.  "Mo"  denotes  "molecule".  "Ra"  denotes  "radi¬ 
cal". 


Fig.  2:  Kc  versus  6  for  po  =  20  bar.  Symbols  □, 
>,  A,  Oj  V,  <1  represent  selected  data  from  the  LLNL  runs; 

lines  •••,—  •  —  •  — , - , - , - — 

represent  the  corresponding  fits.  Temperatures  units  are  in 

degrees  K.  The  legend  is:  •••</>  =  0.5,  Ts  =  825; - 

4>  =  0.5, Ta  =  715,  —  •  —  •  —  </>=  1.0, Ts  =  825; - 

4>=l.0,Ts  =  715;  —  <f>  =  2.0,  Ts  =  825;  - 
4>  =  2.0,  Ts  =  715. 


b  0 

Fig.  3:  Kc  versus  6  for  p  =  1/3  at  different  po  and 
Ts.  Symbols  □,  >,A,0>V,<  represent  selected  data 

from  the  LLNL  runs;  lines  - , - , 

- —  represent  the  corresponding  fits. 

Units:  pressures  in  bar,  temperatures  in  degrees  K.  The  leg¬ 
end  is:  •  •  •  po  =  10,  Ts  =  825; po  =  10,  Ts  =  715, 

-•-•-po  =  20,  Ts  =  825; - po  =  20,  Ts  =  715; 

---po  =  40,  Ts  =  825; -••-••- po  =  40, Ts  =  715. 


Fig.  4:  Kc  versus  6  for  a  variety  of  conditions.  Symbols  □ , 
>,  A,  Oi  V,  <  represent  selected  data  from  the  LLNL  runs; 

lines  •••,  — , - , - ,  —  —  represent 

the  corresponding  fits.  Units:  pressures  in  bar,  temperatures 

in  degrees  K.  (a)  <j>  =  4,  (•  •  •  po  =  10,  Ta  =  825; - 

po  =  10,  Ta  =  715,  —  •  —  •  —  po  =  20,  Ts  =  825; 

- p0  =  20,  Ts  =  715;  -  -  -  p0  =  40,  Ts  =  825; 

—  —  po  =  40,  Ta  =  715).  (b)po  =  10  and <j>  =  1 

=  825; - Ts  =  715,  -  •  -  •  -  Ts  =  655; 

- Ta  =  560). 
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Modeling  of  Alkane  Oxidation  using  Constituents  and 

Species 
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Jet  Propulsion  Laboratory * ,  California  Institute  of  Technology, 

Pasadena  CA  91109-8099 

California  Institute  of  Technology** ,  Mechanical  Engineering ,  Pasadena,  CA  91125 

A  chemical  kinetics  reduction  model  is  proposed  for  alkane  oxidation  in  air  that  is  based 
on  a  parallel  methodology  to  that  used  in  turbulence  modeling  in  the  context  of  Large 
Eddy  Simulation.  The  objective  of  kinetic  modeling  is  to  predict  the  heat  release  and 
temperature  evolution.  In  an  a  priori  step,  a  categorization  of  time  scales  is  first  conducted 
to  identify  scales  that  must  be  modeled  and  scales  that  must  be  computed  using  progress 
variables  based  on  the  model  for  the  other  scales.  First,  a  decomposition  of  heavy  (carbon 
number  greater  or  equal  to  3)  hydrocarbons  into  constituents  is  proposed.  Examination  of 
results  obtained  using  the  LLNL  heptane-oxidation  database  in  conjunction  with  Chemkin 
II  shows  that  (i)  with  appropriate  scaling,  the  total  constituent  mole  fraction  behaves  in 
a  self-similar  manner  and  the  total  constituent  molar  density  rate  follows  a  quasi-steady 
behavior,  and  (ii)  the  light  species  can  be  partitioned  into  two  subsets  according  to  whether 
they  are  quasi-steady  (nine  species)  or  unsteady  (11  species).  The  twelve  progress  variables 
represented  by  the  total  constituent  molar  density  and  the  molar  densities  of  the  unsteady 
light  species  are  defined  to  be  a  base  from  which  the  system’s  behavior  can  be  reproduced. 

This  is  a  dramatic  reduction  from  the  160  species  (progress  variables)  and  1540  reactions 
in  the  LLNL  set  to  12  progress  variables,  16  quasi-steady  rates  (associated  with  heavy 
species),  162  conventional  reaction  rates  (light  species)  and  11  other  functional  forms  (i.e. 
fits  for  the  mean  heavy-species  heat  capacity  at  constant  pressure,  the  enthalpy  release  rate 
of  the  heavy  species,  and  the  molar  fraction  of  quasi-steady  light  species).  A  summary  of 
the  model  is  presented  explaining  the  curve  fits  that  constitute  the  model,  namely  (1)  for 
the  constituent  molar  density  rate  along  with  the  corresponding  enthalpy  production  rate, 

(2)  for  the  quasi-steady  species  mole  fraction,  and  (3)  for  the  contribution  from  the  heavy 
species  to  the  unsteady  light  species  reaction  rates.  The  proposed  kinetic  mechanism  is 
valid  over  a  pressure  range  from  atmospheric  to  60  bar,  temperatures  from  600  K  to  2500 
K  and  equivalence  ratios  from  0.125  to  8.  This  range  encompasses  diesel,  HCCI  and  gas 
turbine  engines,  including  cold  ignition;  and  NOx,  CO  and  soot  pollutant  formation  in  the 
lean  and  rich  regimes,  respectively.  Highlights  of  the  a  priori  model  results  are  illustrated 
for  a  variety  of  initial  conditions.  Results  from  a  posteriori  tests  are  shown  in  which  the 
model  predictions  for  the  unsteady  light  species  and  the  temperature  are  compared  to  the 
equivalent  quantities  based  on  the  LLNL  dataset. 

I.  Introduction 

The  challenge  of  modeling  turbulent  reactive  flows  is  so  considerable  that  the  activity  has  traditionally 
been  decomposed  in  its  two  essential  parts:  kinetics  and  turbulence.  Traditionally,  modeling  of  chemical 
kinetics  has  proceeded  on  a  separate  path  from  that  of  turbulence  which  also  includes  canonical  models  for 
turbulence/reaction  interaction.  The  only  constraint  to  kinetic  modeling  was  that  it  should  be  compact 
enough  to  be  computationally  efficient  when  included  in  a  complex  turbulent  combustion  code.  However, 
there  are  definite  advantages  on  approaching  chemical  kinetic  modeling  in  a  similar  manner  to  turbulent,  flow 
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modeling  because  if  the  concepts  are  similar,  the  hope  is  that  the  models  will  mesh  better  and  the  results 
will  be  easier  to  understand.  This  is  the  approach  taken  here.  The  spirit  of  the  chemical  kinetic  modeling 
approach  is  that  of  Large  Eddy  Simulations  (LES)  in  which  kinematically  energetically-significant  flow  scales 
are  computed  and  the  others  are  modeled;  in  turbulence,  the  large  flow  scales  constitute  the  former  category 
and  the  small  flow  scales  constitute  the  latter  category.  The  chemical  kinetics  parallel  is  to  obtain  a  model  in 
which  one  retains  only  the  thermodynamically  energetically-significant  chemical  scales  as  progress  variables 
and  models  the  fate  of  the  other  scales.  But  the  parallel  approach  between  kinetics  and  turbulence  can 
go  even  further.  In  turbulence  modeling,  a  common  methodology  is  to  assess  the  behavior  of  the  modeled 
scales  by  analyzing  databases  created  using  Direct  Numerical  Simulations  (DNS)  in  which  all  flow  scales  are 
computed;  indeed,  current  experimental  diagnostics  do  not  permit  the  same  thoroughness  of  information  as 
that  obtained  from  DNS.  The  DNS  data  is  analyzed  in  what  is  called  an  a  priori  study  to  inquire  about  the 
behavior  of  the  small  scales  and  propose  mathematical  forms  which  fit  this  behavior.  The  a  priori  study 
is  followed  by  an  a  posteriori  study  where  the  proposed  mathematical  forms  are  inserted  into  the  model  to 
evaluate  its  performance  when  compared  to  the  DNS  database  at  the  LES  resolution. 

A  complete  analogy  between  turbulence  and  kinetics  can  be  made  by  observing  that  kinetic  elemental 
or  skeletal  mechanisms  can  serve  for  reduced  kinetics  the  role  that  DNS  serves  for  LES,  in  which  case 
reduced  kinetic  mechanisms  can  be  viewed  as  the  complement  to  LES  in  achieving  the  goal  of  accurate 
computationally-efEcient  turbulent  reactive  flow  simulations.  The  analogy  between  reduced  kinetic  models 
and  LES  is  not  entirely  surprising  since  each  chemical  species  has  a  characteristic  time  scale  and  in  the 
kinetic  reduction  it  is  desirable  to  compute  only  those  entities  (e.g.  species,  combination  of  species,  radicals, 
combination  of  radicals,  etc.)  having  an  essential  characteristic  time  scale  (to  be  defined)  and  model  the 
kinetics  of  the  remaining  entities. 

Thus,  there  are  two  important  components  to  this  study,  namely  the  a  priori  analysis  and  the  a  posteriori 
evaluation.  Whereas  turbulence  modeling  benefits  from  decades  of  work  using  the  DNS/LES  concept  which 
originated  in  atmospheric  turbulence  predictions  in  the  1960s,  the  present  work  is  the  first  investigation  to 
take  this  approach  in  kinetic  reduction  modeling.  Therefore,  it  is  first  necessary  to  produce  a  categorization 
of  scales  analogous  to  the  large  and  small  scales  of  turbulence,  then  it  is  required  to  propose  mathematical 
forms  for  the  scales  that  will  be  modeled  rather  than  computed  as  progress  variables  in  the  reduced  kinetic 
model,  and  finally  it  is  required  to  perform  an  a  posteriori  study  in  order  to  evaluate  the  chemical  kinetic 
model  versus  the  elemental  or  skeletal  mechanism  for  those  species  predicted  by  the  reduced  model.  The 
present  model  is  for  a  constant-volume  situation,  so  as  to  be  consistent  with  the  requirement  of  a  LES  grid. 

In  Section  II  we  propose  such  a  categorization  through  the  definition  of  a  total  constituent  radical  and 
the  partition  of  the  light  species  set  into  a  set  of  modeled  species  and  a  set  of  progress  variable  species.  That 
is,  rather  than  following  all  species  through  their  reaction  coordinates,  we  follow  a  reduced  set  of  reaction 
coordinates  (i.e.  progress  variables);  this  reduced  set  is  called  a  base.  Section  III  is  devoted  to  the  a  priori 
model.  In  Section  IV,  the  behavior  of  the  base  set  is  briefly  examined  a  priori  in  the  context  of  a  detailed 
reaction  model  for  n-heptane  combustion  as  given  by  LLNL,1  and  the  proposed  mathematical  forms  are 
tested  against  the  LLNL  kinetic  scheme;  some  a  priori  results  were  previously  described  in  more  detail2  and 
we  only  give  here  highlights  of  the  model  and  results  that  were  not  previously  available.  Section  V  is  devoted 
to  a  limited  a  posteriori  evaluation  of  the  model.  Finally,  in  Section  VI  a  summary  of  the  accomplishments 
and  a  perspective  on  future  work  are  offered. 

II.  Categorization  of  chemical  scales 

The  categorization  of  scales  into  modeled  and  computed  is  produced  for  the  primary  objective  of  the 
accurate  determination  of  the  energetics,  i.e.  heat  release  and  temperature  evolution;  at  this  initial  stage 
of  model  development,  neither  NCL  nor  soot  kinetics  are  addressed.  All  species  involved  in  hydrocarbon 
oxidation  are  partitioned  into  light  species  or  heavy  species;  the  heavy  species  are  those  having  a  carbon 
number  n  fs  3.  The  proposed  model  is  for  alkane  combustion  using  as  a  template  the  n-heptane  LLNL1  data 
mentioned  above.  Thorough  examination  of  the  results  using  the  data  led  to  the  following  proposed  repre¬ 
sentation.  The  heavy  species’  rate  evolution  is  modelled,  resulting  in  a  single  progress  variable  constructed 
as  the  total  constituent  radical  molar  density  of  the  heavy  species,  as  defined  in  Section  IIIA.  The  light 
species  set  is  composed  of  two  subsets:  the  first  subset  is  that  of  species  characterized  by  a  quasi-steady 
behavior  which  is  modeled,  and  the  second  subset  is  that  of  unsteady  species  which  are  computed  as  progress 
variables,  both  of  which  are  described  in  Section  IIIB.  That  is,  rather  than  computing  all  species  through 


2  of  14 


American  Institute  of  Aeronautics  and  Astronautics 


their  reaction  coordinates,  a  reduced  set  of  reaction  coordinates  is  computed  that  contains  an  aggregate 
resulting  from  the  heavy  species  decomposition  and  the  unsteady  light  species.  This  reduced  set  is  called  a 
base. 


III.  A  priori  modeling  and  the  resulting  conservation  equations 

The  a  priori  model  addresses  all  heavy  species  and  some  of  the  light  species  (i.e.  all  quasi-steady  light 
species). 


A.  Heavy  species:  the  total  constituent  model  as  progress  variable 

Consider  the  count  of  atoms  (H,C,0,N)  in  species  i.  We  define  Wji  as  the  number  of  species  j  in  species  i 
that  produces  the  correct  atom  count,  where  the  index  j  denotes  the  species  set  of  air  (02,N2)  and  final 
combustion  products  (H20,C02).  (This  means  that  mathematically,  we  consider  species  i  as  being  composed 
of  air  and  final  combustion  product  species.)  The  molar  densities  of  species  in  this  set  are 


•Vi  =  -V"  W 

where  N J  is  the  effective  total  molar  density  of  species  j.  The  same  exercise  of  atom  counting  and  species 
representation  by  the  j  set  can  be  performed  for  all  species  in  the  mixture  to  obtain  Nj.  That  is,  the  atomic 
elements  molar  densities  are  represented  by  species  j  molar  densities.  Conservation  of  atoms  in  the  reactions 
makes  unnecessary  utilizing  rate  expressions  for  reaction  coordinates  of  air  species  and  complete-combustion 
products,  as  a  one-to-one  map  exists  between  these  and  the  atom  set.  The  heavy  species  constituent  radicals 
(CH2,  CH3,  CH,  C2H3,  C2H2,  C2,  HC2,  CO  (keto),  HCO,  HO,  H02,  00,  O)  form  a  set  of  13  entities  from 
which  any  heavy  species  or  radical  may  be  constructed.  The  total  constituent  molar  density  is 

13 

Nc  =  ^Nk  (2) 

k=  1 


where  the  molar  density  of  constituent  k  is  Nk-  The  decomposition  into  constituents  is  similar  to  group 
additivity,3,4  except  that  it  considers  only  first  order  (compositional)  effects  and  only  interactions  with 
adjacent  groups.  The  reaction  rate  associated  with  Nc  is  defined  by 


I<r  = 


d(\nNc) 

dt 


(3) 


We  define  a  reference  density  for  nitrogen  N2  from  the  dry  air  value  at  pressure  prej  =  1  bar  and  Tref= 
298.15  K,  giving  31.5  mol/m3.  Using  Nref,  a  dimensionless  N2  molar  density  is  defined  as 


N*  =  ^ 
Nref 


(4) 


which  acts  as  a  convenient  surrogate  variable  for  the  pressure  p  as  the  N*  value  is  essentially  rate  invariant. 
That  is,  the  partial  pressure  of  N2  is  the  overwhelming  contribution  to  p.  Basically,  the  ratio  of  N*  to  po  is 
nearly  constant.  A  normalized  temperature  variable  9  is  defined  by 


T  —  Ts 

(5) 

Tr(<t>,N*) 

2065 

(6) 

,  1.5  +  1.310 

(!) - 

(7) 

1  + 0.710 +1.102 

where  Ts  is  a  modeling  parameter  representing  the  lowest  temperature  at  which  Kc  >  0.  The  9  definition  is 
such  that  Nc  decreases  during  stoichiometric  combustion  from  its  initial  value  by  three  orders  of  magnitude 
(delving  into  higher  order  of  magnitude  decrease  runs  the  risk  of  encountering  round-off  and  truncation 
errors)  for  9  >  0.6,  a  value  chosen  so  that  all  9  values  remain  below  unity  for  all  test  case  calculations. 
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B.  Light  species:  a  modeled  subset  and  a  progress  variable  subset 

Light  species  are  not  subject  to  meaningful  decomposition  according  to  the  above  procedure.  All  light  species 
are  listed  in  Table  1.  Examination  of  runs  made  using  the  LLNL  database  in  Chemkin  II  indicates  that  the 
light  species  should  be  categorized  in  two  subsets:  one  which  is  modeled  and  one  which  is  computed  subject 
to  the  heavy  species  model. 


1.  Modelled  light  species 

The  species  in  the  first  subset  are  the  radicals  0,  C'H,  CH2,  CH3,  HO,  HCO,  H02,  HC2,  C2H3  which  have 
a  quasi-steady  behavior  and  require  mathematical  fits  of  their  mole  fraction  value  as  a  function  of  the  state 
variables  ( T,N*,q i)  and  modeling  parameters  (here,  only  Ts).  There  are  nine  of  these  species. 


2.  Progress  variable  light  species 


The  species  in  the  second  subset  require  rate  equations.  This  set  consists  of  H20,  C02,  02,  H,  CO,  H2, 
CH4,  H202,  C2H2,  C2H4,  CH20.  The  contribution  of  the  heavy  group  to  the  rate  of  change  of  a  light  species 
molar  density  iVj  with  mole  fraction  Xi  (relative  to  the  light  group)  is  expressed  in  terms  of  quasi-steady 
gain  and  loss  rates  KG-,  and  KL-,  as 


dNj 

dt 


=  Nc{KGi  -  XiKLi ) 

heavies 


(8) 


where  I\Gi  and  KLi  must  be  modeled  consistently  with  the  heavy  species  model.  There  are  eleven  of  these 
species. 


C.  Model  summary  for  species  and  computation  of  energetics 

The  base  set  is  composed  of  a  global  (i.e.  total)  constituent  radical  and  11  light  molecules  or  free  radicals 

of  the  second  light  species  subset.  Thus,  there  are  only  12  progress  variables.  To  compute  Nc  one  must 
model  Kc.  To  compute  the  20  light  species,  one  must  model  X,:  of  the  first  light  species  subset  and  compute 
the  conventional  light  species  reaction  rates  of  the  second  light  species  subset  for  which  models  of  KGi  and 
KLi  are  needed  consistent  with  the  fact  that  the  heavy  species  do  not  appear  explicitly  but  rather  through 
the  total  constituent  Nc. 

Aside  from  the  molar  reaction  rates,  determining  the  temperature  evolution  in  a  reactive  system  requires 
knowledge  of  the  species  molar  enthalpies  and  heat  capacities.  For  species  i,  the  molar  enthalpy  may  be 
expressed  as  hi  =  h®  +  hi(T )  where  h(-  is  the  heat  of  formation,  hi  is  the  sensible  enthalpy  and  T  is  the 

temperature.  The  heat  of  formation  is  taken  at  the  above  reference  conditions,  giving  hi  =  ^  CpgdT 

where  Cvg  is  the  molar  constant-pressure  heat  capacity.  A  heat  of  combustion  for  species  i  is  given  by 

hi  =  -  Y  wPhJ  (Q) 

3 

under  the  assumption  that  water  remains  in  vapor  state.  The  energy  equation  is 


dT 


NcCp,h+  Y  CrdNi  -£  =  ~  E  htTZt+Nc(RuTref)Kh 


iElights 


where 


IZi  = 

Ru  is  the  universal  gas  constant, 


(dNi' 

\  _  m 

V  dt  , 

reac  dt 

iElights 


heavies 


dNi 

dt 


lights 


cp,h  = 


(^-jlCLheavies  C'pjNl') 

Nr  ’ 


and  the  rate  of  enthalpy  change  due  to  heavy  species  l  chemical  kinetic  changes  is  defined  as 

/  \ 

w  )  rTi  1 

KlEheavies 


Kh  =  -  (  5]  hiHi 


B  T  tN 

rej 1  y  c 


(10) 

(11) 

(12) 

(13) 
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Clearly,  since  Kf,.  is  associated  with  the  heavy  species,  it  must  be  modeled. 

For  the  energetics,  values  of  hc  are  fitted  from  the  literature  and  Cv' s  for  the  light  species  are  obtained 
as 


=  ah  +  bh  In 
Ru 


T 


T, 


ref 


(14) 


Other  properties  of  the  light  molecules  are  given  in  Table  l.5  8 


IV.  Highlights  of  a  priori  results 

A.  Examination  of  the  database  using  the  proposed  scaling 

One  of  the  most  important  findings  of  the  categorization  of  scales  identifying  Nc  as  the  important  dependent 
variable  representing  the  heavy  species  is  that  when  Nc  is  scaled  by  (</>  x  N*)  and  the  evolution  is  followed 
as  a  function  of  9  rather  than  using  T,  a  self  similar  behavior  emerges  for  Nc/(<j>  x  N*),  as  shown  in  Fig. 
1  using  the  LLNL  kinetics.  The  similarity  achieved  with  the  two  normalized  variables  is  remarkable,  and 
only  at  a  combination  of  high  po  values  and  small  Ts  values  (typical  of  the  lowest  T  limit  at  which  Kc  rates 
are  finite),  are  small  to  moderate  departures  from  the  excellent  normalization  perceptible.  Even  for  as  rich 
a  mixture  as  <f>  =  2,  similarity  holds,  making  this  model  valid  in  realms  beyond  those  in  advanced  reduced 
schemes9  where  the  upper  limit  of  the  scheme  validity  is  (p  =  1-5,  at  most10  </>  =  2.0  or  exceptionally11  (f>  =  3. 
This  self  similar  behavior  can  be  understood  as  a  reduction  in  dimensionality  of  the  problem  but  should  not 
be  confused  with  the  intrinsic  low-dimensional  manifold12  because  that  concept  was  developed  for  actual 
chemical  species,  whereas  our  findings  are  for  Nc.  The  fact  that  Nc  can  embody  the  evolution  of  all  heavy 
species  is  consistent  with,  and  can  be  traced  to,  the  fact  that  as  T  increases,  the  heavy  species  decompose 
and  no  longer  play  a  role;  instead,  the  products  of  this  decomposition  determine  the  reaction  evolution. 

The  above  findings  were  for  the  LLNL  data.  The  a  priori  study  is  meant  to  propose  models,  that  is, 
mathematical  forms,  for  this  behavior.  This  model  is  briefly  described  next. 

B.  Proposed  model 

Examination  of  Fig.  1  shows  that  by  9  ~  [0.55,  0.6],  Nc  ~  0,  meaning  that  modeling  errors  past  that  9  range 
are  somewhat  irrelevant  to  the  model’s  accuracy  since  the  reaction  is  practically  finished.  Although  it  would 
be  tempting  to  curve  fit  dNJdt  from  Fig.  1  and  adopt  it  as  a  representative  kinetic  rate,  the  fact  is  that  the 
scaling  hides  the  precise  detailed  functional  form  of  Kc.  For  consistency,  Kc  should  be  curve-fitted  similar  to 
Kfr  and  other  reaction  rates.  These  Kc  mathematical  forms  and  an  assessment  of  how  well  they  fit  the  LLNL 
data  were  previously  described.2  According  to  eq.  10,  to  recover  the  value  of  T  in  the  reduction  scheme, 
the  heavy  species  model  should  focus  on  an  appropriate  representation  of  Cpj,  and  Kh.  Plots  (not  shown)  of 
Cp,h  at  various  po  values,  calculated  using  the  LLNL  model  over  a  wide  range  of  <p ,  show  that  these  curves 
exhibited  a  very  modest  variation  over  the  range  of  strong  Nc  decay;  moderate  (i.e.  up  to  50%)  variation 
occurs  only  after  Nc  values  are  very  small  to  negligible.  This  indicates  that  the  T  recovery  is  primarily 
governed  by  Kh  as  far  as  heavy  species  modeling  is  concerned.  Comparison  of  the  model  for  Kh,  which  was 
essentially  similar  to  that  for  Kc  showed  that  it  fitted  remarkably  well  the  results  from  the  LLNL  database 
for  a  variety  of  <p  and  Ts,  including  the  difficult  to  model  very  rich  mixtures.2 

The  model  is  completed  with  the  fitting  of  quasi-steady  gain  rates  KGi  (where  i  stands  for  H20,  C02, 
02,  H,  CO,  H2,  CH4,  H202,  C2H2,  C2H4,  CH20),  of  quasi-steady  loss  rate  KLi  (where  i  stands  for  H20, 
02 ,  H,  H2,  this  rate  being  null  for  the  other  light  species)  and  of  the  quasi-steady  light  species  mole  fraction 
Xi  (O,  CH,  CH2,  CH3,  HO,  HCO,  HOa,  HC2,  C2H3).  Selected  plots  for  KGH.2  and  I<GH  are  presented  in 
Fig.  2;  the  ratio  RLi  =  KLi/KGi  is  illustrated  in  Fig.  3  for  i  =  H2  and  i  =  H;  and  X 'oh  is  depicted  in 
Fig.  4.  The  model  for  these  three  quantities  is  constructed  considering  KGi,  RLi  and  Xi  as  functions  of 
T  or  9,  with  N*,<p  and  Ts  as  being  parameters.  Both  KGi  and  X,  are  split  into  a  low  T,  low-rate  portion 
termed  the  incubation  region  (9  <  0.2),  and  a  high  T,  high-rate  portion  termed  the  fast-rate  region.  The 
incubation  region  variations  exhibit  a  maximum  value  of  KGi  or  X$  followed  by  a  dip  to  a  minimum  from 
which  a  continuous  increase  is  observed  to  the  fast  rate  region.  Conversely,  RLi  first  exhibits  a  minimum, 
and  then  a  maximum.  Curve  fits  to  the  incubation  region  results,  obtained  using  Chemkin  II  with  the  LLNL 
data,  are  generated  utilizing  either  one  of  the  two  following  methods.  In  the  first  method,  sometimes  used 
for  KGi  or  X;,  one  uses  a  cubic  transformation,  e.g.  KGi(T )  to  y{T),  such  that  y  =  —  1  at  the  maximum 
point  and  y  =  1  at  the  minimum  point.  Then  values  of  T  at  fixed  y  values  are  fitted  in  terms  of  parameters 
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( N*,(f>,Ts );  these  values  then  generate  the  continuous  curve  y(T)  by  interpolation  of  the  discrete  values. 
In  the  second  method,  the  T  intervals  before  the  maximum  point,  between  maximum  and  minimum  point 
and  after  minimum  point  are  treated  separately.  Each  of  these  intervals  is  divided  into  equal  T  slices,  and 
logarithm  values  at  slice  boundaries  are  functionally  fitted  in  terms  of  parameters  (N* ,<j),Ts).  The  discrete 
set  of  these  equally  spaced  functional  forms  is  used  to  generate,  by  polynomial  interpolation,  the  continuous 
function.  This  second  method  is  also  used  sometimes  for  the  high-rate  region  fits.  Alternately,  in  this 
region,  for  a  fixed  set  of  (N* the  logarithm  of  the  function  at  all  input  discrete  9  values  is  fitted 
as  a  polynomial  of  9.  The  polynomial  coefficients  are  then  fitted  in  terms  of  (N* ,(f>,Ts).  The  choice  of  the 
particular  method  for  any  i  is  determined  by  the  overall  results  obtained  in  matching  the  input  functions 
provided  by  Chemkin  II  using  the  LLNL  database. 

The  results  in  Figs.  2-4  for  H2  and  H,  as  examples,  show  that  even  at  high  pressure,  the  fits  represent 
the  data  resonably  well.  The  fits  in  Figs.  2a  and  2b  at  Ts  =  610  K  display  an  overestimate  of  the  respective 
KG  in  the  region  9  <  0.1;  the  model  is  more  accurate  as  Ts  is  larger,  e.g.  Figs.  2e  and  2f.  For  the  RL  rates, 
shown  for  H2  and  H  as  examples  in  Fig.  3,  misestimates  also  occur  for  lower  Ts  in  the  9  <  0.1  region  (e.g. 
Figs.  3a  and  3b),  that  are  mitigated  by  the  model  at  larger  Ts.  The  discrepancy  betwenn  LLNL  information 
and  model  for  9  >  [0.55  —  0.6]  is  of  no  concern  since  by  that  station  iVc  ~  0  (see  Fig.  1).  Similarly,  modeling 
discrepancies  for  rich  situations  are  also  of  no  concern  past  the  9  station  at  which  Nc  0. 

An  example  of  light  species  molar  fraction  fits  is  shown  for  OH  in  Fig.  4,  as  OH  is  an  important  radical 
which  is  often  identified  with  the  flame  location.  Over  the  entire  range  of  Ts,  the  model  duplicates  very  well 
the  LLNL  results. 


V.  Results  of  the  a  posteriori  study 

The  developed  model  is  assessed  in  the  a  posteriori  study.  The  assessment  involves  comparing  the  results 
of  the  reduced  kinetic  model  for  T  and  the  unsteady  light  species  H20,  CO2,  O2,  H,  CO,  H2,  CH4,  H2O2, 
C2H2,  C2H4,  CH2O  with  the  predictions  of  the  skeletal  LLNL  mechanism  for  n-heptane.  Since  the  LLNL 
database  results  using  Chemkin  II  does  not  compute  Nc,  an  additional  computation  is  here  made  to  calculate 
this  quantity  as  a  diagnostic  for  the  model,  with  the  understanding  that  it  is  not  a  variable  that  will  be  of 
ultimate  interest  in  model  predictions.  The  comparison  is  made  using  the  developed  model  in  conjunction 
with  Chemkin  II.  The  information  utilized  is  the  set  of  LLNL  reaction  rate  data  for  light  species  interaction 
with  other  light  species,  the  added  rates  of  the  unsteady  light  species  due  to  the  heavy  species  through  the 
fitted  KGi  and  RL-i,  and  quasi-steady  light  molar  densities  through  their  curve-fitted  mole  fractions.  An 
integral  part  of  the  a  posteriori  calculation  is  the  computation  of  Nc  from  the  fits  for  Kc  according  to  eq. 
3  and  of  the  T  evolution  according  to  eq.  10  using  the  fits  for  Kh  defined  in  eq.  13  and  CPth  defined  in  eq. 
12.  Assessment  of  the  reduced  model  involves  comparison  of  its  results  with  the  equivalent  results  obtained 
using  the  full  LLNL  skeletal  mechanism.  This  comparison  is  for  T(t),  Nc(t)  and  Ni(t). 

The  initial  conditions  specify  (j),po  and  T0  which  is  slightly  lower  than  Ts.  Preliminary  reduced  model 
results  indicated  a  strong  sensitivity  to  initial  conditions,  particularly  the  value  of  Ts.  This  led  to  a  further 
examination  of  the  skeletal  LLNL  results  near  the  initial  conditions  of  the  model  calculations.  This  exam¬ 
ination  revealed  that  for  an  initial  T  rise  of  ~  10  K,  the  rate  of  T  increase  is  approximately  linear  with 
T;  i.e.  dT/dt  =  K0(T  —  Tap 0)  where  Tap 0  is  the  temperature  (below  To)  at  which  dT/dt  appears  null,  and 
where  Kq  depends  on  Ts  and  <j>.  Comparison  with  the  LLNL  results  shows  that  Ko  =  ro<fi° A  {Ts/Tref)21  with 
vq  =  2.18  x  10~5  s-1  and  Ts  =  Tq  +  12  within  a  deviation  of  ~  1  K.  The  reduced  model  was  modified  by 
using  this  initial  rate  for  the  region  9  <  10~2,  with  a  smooth  transition  at  the  larger  9  just  past  ICE2  to  the 
rate  based  on  the  energy  equation  involving  the  curve  fit  for  Kh- 

Results  from  the  a  posteriori  evaluation  using  the  modified  reduced  model  are  shown  in  Figs.  5-7.  In 
all  T  evolution  plots  of  Fig.  5  it  can  be  seen  that  the  temperature  increase  past  an  initial,  quiescent  region 
occurs  too  early  compared  to  the  LLNL  results.  This  quiescent  region  is  characterized  by  a  small  T  rise  ~  5 
K  and  nearly  constant  Nc,  with  almost  no  light  species  production.  In  a  second,  active  part  of  the  incubation 
region,  T  increases  to  ~  1000  K,  which  results  from  the  alteration  of  the  structure  of  heavy  radicals  (i.e.  their 
decomposition)  the  changes  of  which  are  represented  by  Kh  through  IZi  of  eq.  13  and  by  the  corresponding 
enthalpies.  The  consequence  of  the  T  increase  is  a  Nc  decrease;  light  species  begin  to  form  during  this  time 
period.  The  time  interval  of  this  second  region,  which  is  also  part  of  the  incubation  period,  decreases  with 
increasing  po-  The  end  of  this  second  region  is  termed  the  ignition  time.  A  successful  model  will  predict 
both  the  ignition  time  and  the  end  of  the  quiescent  portion  of  the  incubation  region.  Fig.  5  shows  that 
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with  the  exception  of  the  lowest  pressure,  po  =  10  bar,  all  the  modified  model  simulations  accurately  recover 
this  ignition  time  and  the  subsequent  temperature  rise,  but  none  of  them  recovers  the  end  of  the  quiescent 
incubation  time  which  corresponds  to  9  ~  10~2.  Following  the  small  T  jump  from  the  quiescent  to  the 
active  incubation  region,  9  becomes  10_1;  thus,  it  is  apparent  that  modeling  accuracy  should  be  improved 
for  9  <  10-2.  With  respect  to  the  ignition  time,  the  reduced  model  also  shows  a  very  sensitive  dependency 
on  Ts,  with  substantial  changes  in  ignition  times  when  Ts  changes  as  little  as  3  K.  The  conclusion  is  that 
although  both  LLNL  results  and  our  model  shows  this  extreme  sensitivity  to  Ts,  the  model  in  the  9  <C  1 
region  does  not  accurately  reproduce  the  LLNL  results;  this  is  because  the  current  model  does  not  address 
in  enough  detail  the  region  9  =  O(10~2),  being  based  on  a  relatively  sparse  selection  of  LLNL  results  so  as 
to  minimize  both  the  fitting  effort  and  the  eventual  computational  cost  of  the  model.  At  this  juncture,  it  is 
clear  that  additional  modeling  effort  must  be  devoted  to  the  quiescent  incubation  region.  The  splitting  of 
the  incubation  region  into  the  quiescent  part  followed  by  a  pre-ignition  T  jump  into  a  region  of  relatively 
slow  T  change  -  i.e.  into  the  active  incubation  region  -  may  be  attributed  to  the  form  of  the  quasi-steady 
rate  functions.  These  rates  typically  feature  peaks  at  9  <  10  1 ,  followed  by  a  significant  dip  to  a  minimum, 
before  an  increase  to  the  fast  rate  region.  This  rate  fall-off  coincides  with  the  slow  T  rise  in  the  active 
incubation  region.  The  conjecture  is  that  this  general  functional  dependency  is  the  essence  of  the  Negative 
Temperature  Coefficient  (NTC)  behavior. 

Considering  the  molar  densities  illustrated  in  Figs.  6  and  7,  the  model  qualitatively  reproduces  all  of 
the  timewise  evolutions,  but  the  quantitative  agreement  is  only  good  to  fair.  The  very  good  quantitative 
agreement  of  Nc  versus  9  of  Fig.  6a  hides  the  quantitative  disagreement  of  Nc  versus  t  of  Fig.  6b.  Although 
it  would  seem  that  this  disagreement  is  due  to  an  inaccurate  Kc  fit  for  9  <  10”2,  the  fact  is  that  the 
discrepancy  in  the  times  of  initial  decrease  of  Nc  (as  shown  in  Fig.  6b)  matches  that  of  times  in  Fig.  5a  at 
the  end  of  the  quiescent  region.  This  explains  the  more  favorable  comparison  in  Fig.  6a  which  masks  this 
discrepancy.  Thus,  the  lack  of  quantitative  Nc(t)  rendition  of  the  LLNL  results  is  traced  to  the  necessity  to 
generate  a  more  accurate  Kh  fit  in  the  9  =  0(  1CA2)  region,  as  already  discussed  in  conjunction  with  the  T(t ) 
variation.  The  need  to  capture  this  delay  at  the  beginning  of  the  active  portion  of  the  incubation  region  is 
also  evident  in  the  evolution  of  light  species  as  a  function  of  t  depicted  in  Fig.  7. 

VI.  Conclusions 

A  computationally  efficient  kinetic  reduction  has  been  proposed  for  alkanes  that  has  been  illustrated  for 
n-heptane  using  the  LLNL  heptane  mechanism.  The  model  is  consistent  with  turbulence  modeling  in  that 
scales  were  here  first  categorized  into  either  those  modeled  or  those  computed  as  progress  variables.  Species 
were  identified  to  be  either  light  or  heavy  (C  number  n  ^  3).  The  heavy  species  were  decomposed  into 
defined  13  constituents  and  their  total  molar  density  was  shown,  through  scrutiny  of  the  LLNL  mechanism, 
to  evolve  in  a  quasi-steady  manner.  The  light  species  were  shown  through  examination  of  the  LLNL  database 
to  behave  either  in  quasi-steady  or  unsteady  manner.  The  modeled  scales  are  the  total  constituent  molar 
density  rate  evolution  and  the  molar  density  of  the  quasi-steady  light  species.  The  progress  variables  are 
the  total  constituent  molar  density  and  the  molar  densities  of  the  unsteady  light  species.  The  unsteady 
equations  for  the  light  species  contain  contributions  of  the  type  gain/loss  rates  from  the  heavy  species  that 
are  modeled  consistent  with  the  developed  mathematical  forms  for  the  total  constituent  molar  density  rate 
evolution  since  examination  of  these  gain/loss  rates  shows  that  they  also  have  a  quasi-steady  behavior  with 
a  functional  form  resembling  that  of  the  constituent  rate.  This  finding  highlights  the  fact  that  the  fitting 
technique  provides  a  methodology  which  can  be  repeatedly  used  to  obtain  an  accurate  representation  of  full 
or  skeletal  kinetic  models. 

An  a  priori  study  previously  initiated  was  here  completed  with  the  development  of  mathematical  func¬ 
tional  forms  for  the  quasi-steady  light  species  and  for  the  gain/loss  rates  for  the  unsteady  light  species.  The 
proposed  reduced  model  features  12  ultimate  progress  variables  describing  heptane  oxidation.  The  reduced 
model  requires  fitting  16  quasi-steady  rate  functions,  11  curves  for  light  quasi-steady  mole  fractions,  and 
Cp^h-  Also  required  are  conventional  rates  between  members  of  the  light  group,  which  are  available  in  the 
literature.  Some  representative  results  of  the  functional  forms  were  shown.  Mathematically,  the  shown 
fits  achieved  substantial  accuracy  over  seven  rate  decades  considering  the  fact  that  they  depend  upon  four 
variables. 

The  focus  was  here  on  the  a  posteriori  study  where  the  findings  of  the  kinetically  reduced  mechanism  were 
compared  to  corresponding  results  obtained  from  the  LLNL  database.  The  comparisons  show  that  although 
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the  model  is  qualitatively  correct  and  reproduces  the  ignition  time  for  the  range  of  conditions  investigated, 
more  fine-tuning  is  necessary  for  the  mathematical  forms  used  in  fitting  the  energetics  of  the  constituents 
in  a  very  small  time  interval  near  the  simulation  initiation  in  order  to  obtain  better  quantitative  agreement 
with  the  LLNL  results. 

Assuming  success  with  the  modified  reduced  model,  the  advantage  of  the  modeling  approach  is  clear: 
Because  this  model  is  based  on  the  Nc  rate  rather  than  on  that  of  individual  heavy  species,  even  if  the 
number  of  species  increases  with  increased  C  number  in  the  alkane  group,  providing  that  the  quasi-steady 
rate  aspect  persists,  then  extension  of  this  model  to  higher  alkanes  should  be  conceptually  straightforward; 
although  it  remains  to  be  seen  if  the  functional  fits  would  remain  valid  or  would  require  reconstruction.  Such 
a  study  is  planned  next  by  ‘blindly’  applying  the  n-heptane  model  to  iso-octane  and  evaluate  its  predictions. 
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Mo/Ra  h°  hc  ah  bh 


H2 

02 

n2 

c 

h2o 

co2 

N 

H 

HO 

HOO 

O 

CO 

HCO 

CH4 

ch3 

ch2 

CH 

C2H3 

hc2 

C2 

NO 

N02 

h2o2 

HC'OH 

C2H2 

c2h4 


0.0 

0.0 

0.0 

717 

-241.5 

-393.5 

473 

218.0 

38 

10.5 
249.2 
-110.5 

43.1 
-74.6 
146 
390 
596 
300 
566 
838* 
90.3 

33.2 
-136 
-109 
228 

52.5 


241.5 
0.0 
0.0 
1111 
0.0 
0.0 
473 
339 
159 
131 

249.2 
283 
558 
802 
902 
1025 
1110 
1449 
1474 
1625 
90.3 

33.2 
106 

526.5 
1257 
1323 


3.282 

3.476 

3.388 

2.50 

3.688 

4.690 
2.50 
2.50 
3.385 
4.150 
2.536 
3.426 
4.154 
3.797 
4.440 
3.973 
3.220 
5.1 
4.434 
4.58 
3.533 

4.691 
5.269 
4.27 
5.368 
5.383 


0.400 

0.5663 

0.469 

0.0 

1.217 

1.390 

0.0 

0.0 

0.3637 

1.307 

0.0 

0.4749 

1.2875 

4.305 

2.249 

1.3015 

0.7136 

3.5 

1.404 

0.0 

0.4508 

1.151 

1.880 

2.546 

2.294 

4.676 


*  Alternate  value  of  832  from  the  CRC  tables. 


Table  1.  Thermodynamic  properties  of  molecules  and  free  radicals.  h°  and  hc  (heats  of  formation  and  combus¬ 
tion,  respectively),  in  kJ/mol;  constants  ah  and  bh  for  molar  heat  capacity  in  the  form  Cp/Ru  =  ah  -\-bhln(T/Tref); 
Tref  =  298.15  K.  "Mo"  denotes  "molecule".  "Ra"  denotes  "radical". 
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Figure  1. 
database. 


Similarity  plots  of  parameter  Nc/(<p  x  N*)  versus  6  at  (a)  po  =  20  bar  and  (b)  <p  =  1  using  the  LLNL 
Ts  is  in  K. 
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Figure  2.  KG  versus  6  for  po  =  40  bar  at  different  Ts  values.  Symbols  □,  >,  represent  selected  results 

using  the  LLNL  database;  lines  represent  the  corresponding  fits.  Temperatures  units  are  in  degrees  K.  The 

legend  is:  0  =  1/4;  >, - 0  =  1/3;  A, - 0  =  1/2;  O? - 0  =  1 ;  —  0  =  2;  <0>,  •••</>  =  4. 

(a,c,e)  for  H2  and  (b,d,f)  for  H. 
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Figure  3.  RL  versus  6  for  po  =40  bar  at  different  Ts  values.  Symbols  □,  >,A,  O,  <,0  represent  selected  results 
using  the  LLNL  database;  lines  represent  the  corresponding  fits.  Temperatures  units  are  in  degrees  K.  The 

legend  is:  0  =  1/4;  >, - <f>  =  1/3;  A, - 0  =  1/2;  O? - 0  =  1;  <],  —  —  0  =  2;  •<),  •••  0  =  4. 

(a,c,e)  for  H2  and  (b,d,f)  for  H. 
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Figure  4.  Xqh  versus  6  for  po  =  40  bar  at  different  Ts  values.  Symbols  □,  >,  A,  Q,  < ,  ^  represent  selected  results 
using  the  LLNL  database;  lines  represent  the  corresponding  fits.  Temperatures  units  are  in  degrees  K.  The 
legend  is:  0  =  1/4;  >, <f>  =  1/3;  A, <j>  =  1/2;  O? - </>  =  1 ;  =  2; -Os  •••</>  =  4. 


t  (s) 


b 


t(s) 


Figure  5.  Timewise  evolution  of  T  at  </>  =  1  and  To  =  820  K  at  several  po  :  (a)  10  bar,  (b)  20  bar,  (c)  40  bar  and 
(d)  60  bar.  The  line  is  the  model  computation  and  the  symbols  are  obtained  using  the  LLNL  database. 
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Figure  6.  Evolution  of  Nc  as  a  function  of  (a)  6  and  (b)  t  for  (j)  =  1,  po  =  10  bar  and  To  =  820  K.  The  lines 
correspond  to  the  model  predictions  and  the  symbols  are  obtained  using  the  LLNL  data. 


Figure  7.  Timewise  evolution  of  molar  densities  of  several  light  species  at  (j>  =  l,po  =  10  bar  and  To  =  820  K. 
(a)  CH4,  (b)  C2H2  and  (c)  CO.  The  lines  correspond  to  the  model  predictions  and  the  symbols  are  obtained 
using  the  LLNL  data. 
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